CYGNSS Surface Heat Flux Product Development

Abstract

Ocean surface heat fluxes play a significant role in the genesis and evolution of various marine-based atmospheric phenomena, from the synoptic scale down to the microscale. While in-situ measurements from buoys and flux towers will continue to be the standard in regard to surface heat flux estimates, they commonly have significant gaps in temporal and spatial coverage. Previous and current satellite missions have filled these gaps; though they may not observe the fluxes directly, they can measure the variables needed (wind speed, temperature and humidity) to estimate latent and sensible heat fluxes. However, current remote sensing instruments have their own limitations, such as infrequent coverage, signals attenuated by precipitation or both. The Cyclone Global Navigation Satellite System (CYGNSS) mission overcomes these limitations over the tropical and subtropical oceans by providing improved coverage in nearly all weather conditions. While CYGNSS (Level 2) primarily estimates surface winds, when coupled with observations or estimates of temperature and humidity from reanalysis data, it can provide estimates of latent and sensible heat fluxes along its orbit. This paper describes the development of the Surface Heat Flux Product for the CYGNSS mission, its current results and expected improvements and changes in future releases.

1. Introduction

Latent and sensible heat fluxes (LHF and SHF, respectively) over the Earth’s oceans, produced by turbulent flow within the planetary boundary layer, can have a significant effect on the genesis and evolution of various weather and climate systems. LHF and SHF are primarily driven by surface wind speeds and differences in temperature and specific humidity between the lowest levels of the atmosphere (~10 meters above the surface) and the sea surface [1]. Influxes of moisture and thermal energy can decrease the static stability within the boundary layer, decreasing the height of the lifted condensation level, potentially leading to the development of dry and moist convection [2,3]. The energy transported between the ocean surface and lower atmosphere influences the development of various weather systems from the microscale up to the synoptic scale, including—isolated and organized convection [4], tropical [5] and extratropical cyclones [6] and large-scale waves (e.g. Madden-Julian Oscillation (MJO)) [7].

There are several methods that can be used to estimate surface heat fluxes but the most common method utilizes the bulk aerodynamic formulas. These formulas relate the turbulent fluxes to the observable spatial and temporal averages [1,8] and can be written as follows:

(1)$LHF= {\rho }_a{L}_v{C}_{DE}U\left(q_s-q_a\right),$

(2)$SHF= {\rho }_a{c}_p{C}_{DH}U\left(T_s-T_a\right),$

Here, ρ_a is the air density at the surface [kg m⁻³]; L_v is the latent heat of condensation (2.5 × 10⁶ J kg⁻¹); and c_p is specific heat at constant pressure (1004 J K⁻¹ kg⁻¹). C_DE and C_DH are, respectively, the exchange coefficients of moisture and sensible heat [unitless]; U is the magnitude of the surface wind speed [m s⁻¹], T_s and q_s are temperature [K] and specific humidity [kg kg⁻¹], respectively, at the surface, while T_a and q_a are the same but at 10 meters above the surface.

Given the impacts LHF and SHF have on weather systems at multiple scales and their rapid spatial and temporal variations, routine and widespread measurements of the surface heat fluxes are necessary. In-situ observations from flux towers and buoys have been and will continue to be the standard for LHF and SHF measurements over the world’s oceans, but they are limited spatially and temporally. While spaceborne instruments are currently unable to directly measure surface heat fluxes, they can be used to estimate the components, like temperature, humidity and wind speed, which are needed to estimate LHF and SHF through the bulk aerodynamic formulas [9]. However, these instruments, such as microwave radiometers and scatterometers in a polar-orbit, have their own limitations. Their orbits can cause them to miss large spatial areas and/or feature infrequent sampling, especially over the lower latitudes. Additionally, their measurements are often affected by the presence of precipitation, leading to inaccurate or missing surface measurements. By utilizing its surface wind speed observations (coupled with estimates of thermodynamic variables from other sources), the Cyclone Global Navigation Satellite System (CYGNSS) can provide improved coverage of surface heat fluxes over the tropical and subtropical oceans with its high temporal and spatial sampling of ocean surface wind speed. Though CYGNSS is in a tropical orbit (35° orbit inclination), it can still observe a large swath of surface heat fluxes over the world’s tropical and subtropical oceans.

The largest latent and sensible heat fluxes are often observed in the extratropical regions during the respective winter seasons of both hemispheres [10]. This is because surface wind speeds and differences in temperature and humidity are consistently at their greatest during the winter. These large heat fluxes are typically concentrated near coastlines where warm ocean waters interact with cold and dry air masses originating over continental landmasses, such as the Western Pacific and Western Atlantic Oceans, which can affect the formation of rapidly developing extratropical cyclones (ETCs) [11,12]. While these fluxes are at their maximum at middle and high latitudes, CYGNSS can make consistent observations up to 38° latitude in both hemispheres [13,14]. As Figure 2 from Yu and Weller [10] indicates, some of the highest fluxes typically observed over the world’s oceans occur within CYGNSS’s observational range.

This paper describes the development of a Level-2 (L2) ocean Surface Heat Flux product for the entire CYGNSS mission, which utilizes its L2 wind speed retrievals, coupled with a reanalysis dataset for the thermodynamic variables (temperature and humidity). The goal of this product is to, along with in-situ and other remote sensing observations, aid the scientific community’s understanding of how latent and sensible heat fluxes impact the genesis and evolution of various weather and climate patterns across the globe by utilizing CYGNSS’s frequent sampling over the tropical and subtropical oceans.

The structure of the paper is as follows—Section 2 introduces in detail the data that are used to estimate the surface heat fluxes for the CYGNSS mission. Section 3 discusses the algorithm being used for the product. Section 4 discusses initial results from the CYGNSS Surface Heat Flux Product across the globe and for two case studies. Section 5 focuses on the comparisons with LHF and SHF estimates from buoy data, as well as how differences in the inputs (wind, temperature, humidity) correlated with differences in the fluxes. Section 6 includes a simple uncertainty and sensitivity analysis of the CYGNSS Flux product. Section 7 and Section 8 contain the discussion, summary, conclusions and future uses and development of the product.

2. Data

2.1. Cyclone Global Navigation Satellite System (CYGNSS)

The Cyclone Global Navigation Satellite System (CYGNSS) was launched on 15 December 2016; it consists of a constellation of eight small satellites in an orbital inclination of 35° that are designed to estimate surface winds over the tropical and subtropical oceans by utilizing Global Navigation Satellite System-Reflectometry (GNSS-R). This technology measures the direct signal from the existing Global Positioning System (GPS) satellite constellation through a zenith antenna and the reflected GPS signal from the ocean surface through its two downward pointing antennae [13].

Every CYGNSS observatory is able to observe up to four simultaneous specular points per second (32 specular points per second for the entire constellation). The scattering map from the reflected GPS signal around the specular point has coordinates of code chip delay and Doppler shift and as such is referred to as a Delay Doppler Map (DDM). In order to estimate the surface wind speeds, the average reflected power around the specular point (DDM average; DDMA) and the slope of the DDM waveform in timed delay coordinates (Leading Edge Slope; LES) are used separately. The DDMA and LES winds are then optimally combined using a Minimum Variance (MV) estimator to produce a single best-estimate wind speed product [13,15]. Time averaging is applied to consecutive samples (DDMs) in order to produce a consistent 25 km spatial resolution for each data product; the number of samples used depends on the incidence angle of the specular point.

This wind speed estimate assumes that the sea state is in equilibrium with the wind speed, formally referred to as “fully developed seas” (FDS). Though these wind estimates are reliable most of the time throughout the tropical and subtropical oceans, they can often be inaccurate at higher wind speeds and in rapidly developing storm systems. Therefore, an alternative wind speed product, “young seas with limited fetch” (YSLF) is produced and designed to capture the sea state when it is not in equilibrium with the local wind. In the YSLF product, the CYGNSS LES and DDMA data are matched with observations from the Stepped Frequency Microwave Radiometer (SFMR) onboard the NOAA P-3 Hurricane Hunters at high wind speeds, typically within tropical cyclones. Under YSLF, DDMA-based wind alone is used, as it has greater sensitivity at higher wind speeds. Additionally, unlike the FDS wind speeds, time averaging is not applied to consecutive samples in order to retain the highest possible horizontal spatial resolution in high wind situations [16]. As one will see in Section 3, however, this does result in noisier data compared to the FDS wind speed estimates. However, the temporal and spatial resolution remain the same between the FDS and YSLF wind speed products from CYGNSS.

The CYGNSS L2 retrieval algorithm produces both the FDS and YSLF wind speed products for the entire mission. Since it does not switch between the products when transitioning between steady winds and winds in rapidly developing systems, it is therefore up to the user to decide when it is best to use the FDS winds and when to use the YSLF winds. To be consistent with the standard CYGNSS L2 products, the inaugural version of the CYGNSS Surface Heat Flux Product also estimates LHF and SHF using both the FDS wind speeds and the YSLF wind speeds separately. It is suggested that the FDS version of wind speeds and fluxes be used across most of the globe and to only use the YSLF versions for high wind and rapidly developing systems where the ocean state most likely does not match the wind speed.

2.2. MERRA-2 Reanalysis Data

The bulk flux algorithm (which will be discussed in the next section) requires inputs of sea surface and near surface thermodynamic variables. These variables, as seen in Equations (1) and (2), include sea surface temperature (T_s), 10-meter air temperature (T_a) and 10-meter specific humidity (q_a). Within the algorithm, surface air pressure (p) is used to estimate air density (ρ) and specific humidity (q_s) is estimated as a function to T_s and p. Given the CYGNSS only provides the wind speeds (U), which is also used to estimate the drag coefficients (C_D), another source is needed for these variables. For the initial version of the product, the values for these variables are obtained from the NASA Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) [17]. MERRA-2 uses data assimilation to combine all available in-situ and satellite observation data with an estimate of the atmospheric state, provided by a global atmospheric model. While there are other products available, MERRA-2 is a good source for this initial version of the Surface Heat Flux Product as it provides reliable data with relatively high spatial and temporal resolution for all the variables needed to estimate LHF and SHF. MERRA-2 offers a temporal resolution of 1-hour and a spatial resolution of 0.5° × 0.625°. Since its spatial resolution is similar to that of the CYGNSS L2 wind speed data (25 km ~ 0.25°), we are able to use a simple nearest neighbor approach in time and space in order to match the CYGNSS specular points to the nearest MERRA-2 grid point.

3. Algorithm Description

3.1. COARE 3.5 Background

Version 3.5 of the Coupled Ocean-Atmosphere Response Experiment (COARE 3.5) is a widely used parameterization of the latent and sensible heat flux transfer coefficients [18]. COARE was initially designed to analyze and understand the processes that occur in the region of the tropical western Pacific Ocean warm pool, such as ocean-atmosphere coupling and multi-scale interactions that extend oceanic and atmospheric influences into other regions [19]. The COARE algorithm is based on Monin-Obukhov similarity theory (MOST) and has been widely used to estimate surface heat fluxes over the oceans [20,21]. The COARE algorithm was initially designed to produce flux estimates in the presence of low to moderate wind speeds, but the latest version has been validated versus direct in-situ flux measurements for wind speeds up to 25 m/s [18].

COARE 3.5 parameterizes the drag coefficients present in the bulk aerodynamic formulas, C_DE and C_DH, as a function of gustiness, surface roughness and atmospheric stability. The parameterization has been derived from direct covariance flux measurements collected from various field campaigns over the midlatitude oceans [18]. The drag coefficient, C_D, within COARE 3.5 is expressed mathematically as:

(3)$C_D\left(\frac{z}{z_o},\frac{z}{L},G\right)=\frac{-\overline{uw}}{U_r{S}_r}=\frac{-\overline{uw}}{U_r^{2}G}={\left[\frac{\kappa }{\ln \left(\frac{z}{z_o}\right)-{\psi }_m\left(\frac{z}{L}\right)}\right]}^2,$

Here, z is the height above the surface [m]; κ is the von Kármán constant (set to a value of 0.4; unitless); z_o is the aerodynamic roughness length [m] and ψ_m is a dimensionless function that accounts for the effects of atmospheric stratification. G is the gustiness parameter, defined as the ratio of the wind speed, S_r [m s⁻¹], to the vector-averaged wind, U_r [m s⁻¹] and uw are the fluctuating horizontal and vertical velocity components [m s⁻¹], respectively, with the overbar representing a time averaging of these components [18,22]. This attempts to account for mass, momentum and heat transfer at low wind speeds, which is non-zero because of gustiness; wind gusts produce shear-driven turbulence that can drive a significant portion of the surface-atmosphere exchanges in convective conditions [18,20].

As mentioned earlier, COARE 3.5 is validated for wind speeds up to 25 m/s. When wind speeds exceed this limit, sea spray ejected from the ocean surface under high-wind conditions has a non-negligible effect on the air-sea fluxes that are observed [23], which is not accounted for within the COARE algorithm. Additionally, the estimation of the drag coefficient begins to break down when wind speeds exceed 25 m/s; these result in increasing errors in the LHF and SHF estimates derived through the bulk aerodynamic formulas. Given these limitations, latent and sensible heat fluxes that are estimated when wind speeds exceed 25 m/s are flagged within the product.

3.2. Algorithm Flow

As mentioned in Section 2.2, the CYGNSS specular points are matched with a MERRA-2 grid point using a simple nearest-neighbor method, as MERRA-2 provides reasonably high temporal (hourly) and spatial (0.5° × 0.625°) resolution output. From here, the matched CYGNSS specular points and MERRA-2 data (MERRA-2 input variables listed in Section 2.2) are ingested into the COARE 3.5 algorithm. As mentioned in Section 2.1, two CYGNSS wind speed products are used—the minimum error variance estimate of the FDS wind speed and YSLF wind speeds. Both wind speed products are used separately in order to calculate two LHF and two SHF products, allowing users to maintain consistency with analyses of CYGNSS surface fluxes and wind speed observations. As in their L2 wind speed counterparts, LHF and SHF estimates computed using FDS winds are recommended for most global measurements, while LHF and SHF estimates computed using YSLF winds are recommended for rapidly developing systems and regions with strong curvature in the flow (i.e. tropical and extratropical cyclones). Currently, CYGNSS does not provide a mixed FDS and YSLF product that switches between the two automatically; therefore, a mixed FDS and YSLF product is not provided for the CYGNSS Fluxes. This is something we hope to address in future versions of the CYGNSS Surface Heat Flux Product.

In the initial data releases of the CYGNSS L2 products (versions 2.1 and previous), there were uncertainties in the Block II-F GPS transmitter power and antenna gain patterns, which significantly and negatively affected the wind speed products. As a result, specular points that use the Block II-F satellites are discarded in the L2 wind speed product [24] and are subsequently flagged in the surface heat flux product. Additional quality flags are produced for LHF and SHF estimates when the range corrected gain (RCG) is less than three, since wind speed estimates can be unreliable at RCG < 3 [25]. This often includes specular points poleward of the 40^th parallel in both hemispheres, as they often have high incidence angles and noisier results, coinciding with low RCG values. Finally, LHF and SHF results are flagged whenever the FDS or YSLF winds exceed 25 m/s for the reasons mentioned in Section 3.1.

Once the inputs from MERRA-2 and CYGNSS are inserted into the COARE algorithm, it produces a first guess of the surface heat fluxes in order to initialize a stability iteration loop. Within this loop, COARE computes the Monin-Obukhov length, roughness length and transfer coefficients along with a stability dependence. For the initial version of the Surface Heat Flux Product, this loop is repeated at least 10 times. This is based on various trials that have shown that this is the minimum needed for the values to reach an asymptote at each specular point. However, future versions of the CYGNSS Fluxes will include more stringent criteria for when this iteration loop should terminate at each specular point. These transfer coefficients from the stability interaction loop are then used to calculate the latent and sensible heat fluxes using the bulk aerodynamic formulas (Equations (1) and (2)) and are repeated for both wind speed products from CYGNSS, as discussed in Section 2.1. For more information regarding the COARE algorithm and its details, readers are encouraged to see Edson et al. [18].

Two versions of latent and sensible heat fluxes are produced within this algorithm, one that uses CYGNSS FDS winds and one with YSLF winds, outputting four total products from this algorithm. These results are flagged if they meet any of the criteria listed earlier in this section. Along with these, an uncertainty parameter is also produced that correlates with the uncertainty analysis that is discussed in Section 6. Finally, variables from MERRA-2, as well as a sample index that matches that of the L2 CYGNSS wind product, is produced so that users can reproduce these fluxes if they so desire.

4. Results

Figure 1; Figure 2 depict two separate full days (split into 12-hour increments) of latent and sensible heat flux estimates from CYGNSS using the FDS wind speeds. On 1 January 2018 (Figure 1), one can see large latent and sensible heat fluxes, primarily along the Western Pacific and Western Atlantic Oceans in the Northern Hemisphere. This is expected, as large air-sea temperature and humidity differences observed during the winter seasons lead to higher heat fluxes, as well as high wind scenarios associated with possible developing winter storms. Meanwhile, the distribution of fluxes are reversed in the 1 July 2018 case (Figure 2) as it is the winter season in the southern hemisphere, leading to higher surface heat fluxes in the South Pacific, Indian and Atlantic Oceans. Strong gradients in mean sea level pressure are co-located with some regions of high heat flux, indicating the correlation between high wind speeds and large positive surface heat fluxes.

In September 2018, Hurricane Florence made landfall along the North Carolina shore but after landfall, its center of circulation moved to the southwest along the coast, causing heavy rainfall and major coastal flooding in the area. As post-landfall Florence travelled along the North Carolina coast, CYGNSS was able to observe the cyclone and estimate its associated surface heat fluxes. Around 18:00 UTC on 14 September, LHF observations from both FDS and YSLF products exceeded 300 W/m² (Figure 3a,b), though SHF observations hardly exceeded 30 W/m² (Figure 3c,d). While there were latent heat flux values in and around Florence up to 300 W/m² prior to landfall, regions containing high fluxes were not as widespread as high-flux regions observed after landfall. While the high wind speeds observed in Florence certainly contributed to the high latent heat fluxes, the large air-sea differences after landfall had an influence as well. In its western half, Florence was able to pull in drier air from the continent, which then interacted with the warm Gulf Stream ocean surface on its southern and eastern halves, contributing to the high latent heat fluxes observed in concert with its higher wind speeds. However, SHF remained low, indicating that while the humidity differences were large, the same could not be said about the temperature differences between the surface and 10 meters. As expected, the surface fluxes that were estimated utilizing the YSLF winds were relatively high; however, one should note that there are fewer points compared to the fluxes estimated with FDS winds. This can be noticed in the southern half of the cyclone. This is expected, as YSLF wind speeds are generally higher than FDS wind speeds and are therefore more likely to surpass the 25 m/s speed limit of the COARE 3.5 algorithm and be removed through the quality flag controls.

Given CYGNSS’s ability to observe low-latitude extratropical cyclones [14], we can also obtain LHF and SHF estimates in extratropical cyclones (ETCs) that develop within CYGNSS’s range. One of the ETCs that CYGNSS observed was a rapidly developing ‘bomb cyclone’ (an ETC in which the mean sea level pressure decreases by at least 24 hPa within 24 hours [26]) along the East Coast of the United States in January 2018 (Figure 4). Though the ETC developed and moved poleward rapidly, CYGNSS was able to observe the equatorward side of the cyclone. The highest surface heat fluxes are typically observed on this side of a marine-based ETC, as strong winds pull cold and dry air from the continent over the warmer and moister ocean surface [6]. As can be seen in this case study on 4 January 2018 around 15z, CYGNSS observed latent heat fluxes exceeding 600 W/m² and sensible heat fluxes of over 300 W/m² as the cyclone continued to strengthen and proceed poleward to the northeast United States. Much like the Hurricane Florence case, the LHF and SHF values estimated with YSLF winds were often higher than those estimated with FDS winds. In this case, one can clearly see that fluxes estimated from the YSLF winds are noisier than LHF and SHF estimated with FDS wind speeds. This is expected because, as mentioned in Section 2.1, time averaging is not applied to consecutive sample points for the YSLF winds, as is done for the FDS winds, so that the highest possible horizontal spatial resolution is obtained for high wind situations.

High surface heat fluxes were observed in both the Hurricane Florence and ‘Bomb Cyclone’ ETC cases, indicating that there is a large amount of energy being transferred from the ocean surface into the lower atmosphere, which could affect how these systems evolved. While previous research has shown that high surface heat fluxes associated with these and other systems can influence their development [4,5,6,7], we cannot speculate at this time exactly how much the fluxes observed in these case studies affected their development solely based on these observations. Future modelling studies could address these specific case studies and examine how the high surface heat fluxes that were observed by CYGNSS affected their evolution.

5. Comparisons with Ground Truth Buoy Data

Direct in-situ measurements of LHF and SHF have been limited during CYGNSS’s mission, as they are often collected only during a small number of field campaigns (e.g. PISTON [27], CAMP2Ex [28]) and on research buoys. While comparisons with additional in-situ data may be obtained for future versions of the CYGNSS Flux Product, at this time we have used surface flux estimates derived from buoy data that are a spatial and temporal match with CYGNSS specular points. Though these buoys do not measure LHF and SHF directly, they do measure wind speeds, temperature and humidity, which can be directly inserted into the same COARE algorithm utilized for the CYGNSS Surface Heat Flux Product.

From 18 March 2017 to 29 January 2019, the CYGNSS Surface Heat Flux Product was validated against the following buoys and buoy networks (Figure 5):

Kuroshio Extension Observatory (KEO) [29]

National Data Buoy Center (NDBC) [30]

Ocean Sustained Interdisciplinary Timeseries Environment observation System (OceanSITES) [31]

Prediction and Research Moored Array in the Tropical Atlantic (PIRATA) [32]

Research Moored Array for African-Asian-Australian Monsoon Analysis and Prediction (RAMA) [33]

Tropical Atmosphere Ocean Array (TAO) [34]

Data from these buoys can be accessed publicly through the National Oceanic and Atmospheric Administration (NOAA) [35,36,37]. Data from the buoys were compared to the CYGNSS Fluxes following Ruf et al. [38]. In the present analysis, CYGNSS derived flux measurements that occurred within 50 km and 0.5 hours of a buoy’s location and observation time were used, which were then collocated by an inverse-weighting scheme [39]. Overall, 83 buoys were used for these matchups, with an aggregated sample size of 21,679 matchups for the FDS products and 20,947 matchups for the YSLF products (after applying CYGNSS Flux quality flags).

5.1. Comparisons with Buoy Data

Figure 6 shows a 2D density scatter plot of the collocated flux samples. The scatter plot demonstrates a good agreement between the CYGNSS Fluxes and the derived buoy fluxes, with a correlation around 0.8 and 0.79 for LHF with FDS and YSLF winds (respectively) and 0.85 for both SHF products (Table 1). The highest density of matchups occur along the diagonal 1:1 line for all CYGNSS Flux products; the highest density of LHF matchups occur between 50-150 W/m², while the highest density of SHF matchups occur around 0 W/m². The CYGNSS Fluxes overall agree well with the estimates from the buoys, though there is a discrepancy at lower values (as seen in Figure 7). However, as the surface heat fluxes increase, there is greater scatter and disagreement between CYGNSS and the buoy data. While some fluxes are overestimated, a majority of the higher CYGNSS Flux values are lower than those estimated from the buoy data. This is likely due to CYGNSS’s increased uncertainty at higher wind speeds, which have yielded lower wind speed values compared to other sources [40]. Since wind speeds correlate with LHF and SHF based on Equations (1) and (2), if CYGNSS winds are underestimating higher wind speed values, then it should also lead to an underestimate of higher latent and sensible heat flux values. However, as will be discussed in the next section, the inputs from MERRA-2 could also be impact the disagreements between CYGNSS and buoy flux measurements.

While the density plots in Figure 6 show that a majority of the fluxes match-up well, there are still significant spread in the data at the lower flux values. Figure 7 is a histogram of the CYGNSS Fluxes and buoy data binned in intervals of 25 W/m². While there are few differences in the SHF histogram, there are larger differences in the LHF plots. For LHF values between 50–200 W/m², we see a greater number of CYGNSS observations compared to the buoy data, while the reverse is true for values less than 50 W/m². These plots help us understand the additional spread around the 1:1 line observed in the density plots (Figure 6). The differences seen in Figure 6 and the scatter observed at higher flux values could be the result of uncertainties in the CYGNSS wind speed observations, as there have been known errors at higher wind speeds [38,40]. However, CYGNSS is likely not the only source of uncertainty and errors at these higher flux values. As Section 5.2 and Section 6 will discuss, differences in variables between MERRA-2 and the buoy data could also be contributing to these differences.

As highlighted by Figure 7, the number of buoy and CYGNSS comparisons at higher flux values were limited in comparison to those at lower LHF and SHF values, which may have resulted in the larger scatter at higher flux values. We expect that future versions of the CYGNSS Fluxes will include comparisons and validation with flux measurements from various field campaigns, as well as a larger set of buoy comparisons, especially at the high fluxes values. In addition to validation with additional sources of data, future releases of the CYGNSS wind speeds are expected to reduce some of the errors at higher wind speeds that are currently present in the current version.

5.2. Flux Component Influence

Figure 8 shows how differences between CYGNSS Flux observations and buoy flux observations change as wind speed measurements from the buoys increases. Here, a positive difference means CYGNSS observations were greater than the buoy data and vice versa for negative differences. At lower wind speed observations, while there is a large spread between the two (similar to Figure 6), based on the quadratic fit the differences between the CYGNSS and buoy fluxes averaged out to be near zero. Meanwhile, as wind speeds from the buoy data increases, the differences between CYGNSS and buoy flux measurements increase, with all four CYGNSS Flux products underestimating LHF and SHF with respect to the buoy data. These differences, as they correlate to wind speed, could explain why we see CYGNSS underestimate the fluxes compared to buoy data in Figure 5. While differences in wind speed could be the main factor in these differences, as CYGNSS winds have been underestimating winds at higher wind speeds [40], we need to look into how all inputs into the flux product could be impacting these results.

Figure 9 shows the flux residual dependency on available bulk flux component (temperature, humidity, wind speed) residual. Here, all differences are calculated as CYGNSS minus buoy and the results are normalized to the range of −10 to 10, similar to previous studies [41]. For both latent heat flux products, differences in wind speed (ΔU) correlate well with ΔLHF. However, though smaller, the Δq_a and Δq_s have a non-negligible impact and have a correlation with ΔLHF as well, indicating that the variables in the MERRA-2 contribute to these differences.

Meanwhile for both products of the sensible heat flux, discrepancies in wind speed between CYGNSS and buoys do not contribute much to the observed ΔSHF. However, ΔT_a and ΔT_s correlate more to the observed ΔSHF. This shows that while discrepancies in CYGNSS and MERRA-2 can equally contribute to observed ΔLHF, it seems that differences in MERRA-2 observations are the main contributor to observed ΔSHF.

One curious oddity in these plots are the large spread in ΔLHF and ΔSHF in the surface humidity and temperature plots, respectively, around Δq_s= −5 and ΔT_s= −5 (originally ΔT_s= 0°C & Δq_s= 0 kg/kg, respectively). While future analysis is needed to address why there is a spread in ΔLHF and ΔSHF when ΔT_s and Δq_s are zero, one hypothesis that could explain these differences are other factors, such as wind and 10-meter humidity/temperature could be playing a much larger role in ΔLHF and ΔSHF at these points. Additionally, it could also be due to the air-sea state or stability effect. However, future analyses are needed to fully address the oddity in this product.

6. Simple Sensitivity and Uncertainty Analysis

While the use of the COARE algorithm for the CYGNSS Flux Product has been effective in estimating LHF and SHF throughout CYGNSS’s mission, inaccuracies and uncertainties in the inputs can result in greater uncertainties in LHF and SHF. As discussed in Section 5, variations in the inputs from CYGNSS and MERRA-2 could significantly alter the results. In order to estimate the uncertainty of the fluxes from these inputs, we first examine the impact of the known and estimated uncertainties of each input variable (wind, temperature, humidity) on the LHF and SHF results.

As an example, an uncertainty value is reported at every specular point for CYGNSS’s FDS and YSLF wind speed product. For this uncertainty analysis, at each specular point, 100 random numbers are created with a mean of zero and a standard deviation equal to the reported wind speed (FDS) uncertainty. These 100 values are then added to the FDS wind speed to produce 100 perturbed wind speed values at that specular point. The same algorithm is then applied to estimate LHF and SHF for each perturbed value of wind speed, keeping the other inputs the same. This produces 100 values each for LHF and SHF at every specular point. Then, the standard deviation of these 100 values are computed to estimate an uncertainty in LHF and SHF associated with the wind speed uncertainty.

This process is repeated for YSLF winds, surface temperature, 10-meter air temperature and relative humidity at each specular point. Unlike CYGNSS, MERRA-2 does not provide uncertainties at each grid point, so an averaged uncertainty for each variable is used across the globe (1°C for air temperature, 0.5°C for SST and 5% for relative humidity [42,43]). Since uncertainties for specific humidity from MERRA-2 were not available, 10-meter specific humidity was converted to relative humidity given its availability [42]. These values are an approximate from similar studies, but future uncertainty analysis and updates for the CYGNSS Fluxes will examine how these uncertainties can vary in different conditions (e.g. precipitation, rapidly developing systems, etc.), as well as include uncertainties in specific humidity at the surface.

Once the standard deviations of LHF and SHF from each variable are computed, they are combined into one standard deviation for LHF and SHF at every specular point using the formula:

(4)${\sigma }_{LHF}=\sqrt{{\sigma }_{LHF,Ta}^2+{\sigma }_{LHF,Ts}^2+{\sigma }_{LHF,RH}^2+{\sigma }_{LHF,FDSwind}^2},$

where σ_LHF,Ta is the standard deviation of LHF when 10-meter air temperature is perturbed, σ_LHF,Ts is the standard deviation of LHF when surface temperature is perturbed, σ_LHF,RH is the standard deviation of LHF when relative humidity is perturbed and σ_LHF,FDSwind is the standard deviation of LHF when FDS wind speed is perturbed (all units of W m⁻²). This same method can be used to estimate σ_SHF, along with σ_{LHF_YSLF} and σ_{SHF_YSLF}, by replacing the last input in Equation (4) with σ_{LHF, YSLFwind} and σ_{SHF, YSLFwind}, respectively.

Given that a single day of CYGNSS observations can yield over 2.5 million specular points throughout the entire constellation, this simple uncertainty analysis was limited to just three days (14-16 September 2017) for this study due to computational and time limits. Future analysis will include more data, as well as analyze how these uncertainties can change in different seasons, environments, basins, and so forth. Specular points were that flagged, under the criteria described in Section 3.2, were not used for the analysis.

As shown in Table 2, the uncertainties at each specular point over this three-day period are averaged to give the overall uncertainty for LHF and SHF, as well as the sensitivity from the individual components.

These values show how the uncertainties from CYGNSS and MERRA-2 can contribute to uncertainties observed in SHF and LHF (11 W m² for both SHF products, 44/51 W m² for the FDS/YSLF versions of LHF product, respectively). As expected, wind speed is the only input whose uncertainty impact changes between the FDS and YSLF products. While there are very small changes for temperature and humidity between the FDS and YSLF products, these are insignificant and can most likely be attributed to differences in sample points (as more YSLF points are filtered out due to higher wind speeds). Despite the large uncertainties associated with the CYGNSS winds, in most products it was not the largest source of uncertainty. Though a 1°C uncertainty was used for 10m air temperature, it led to the largest source of uncertainty in three of the four CYGNSS Flux products. As was seen in Figure 9, differences in CYGNSS and buoy wind speeds did not contribute much to ΔSHF, whereas differences in temperature between MERRA-2 and the buoy data did. This correlates well to the uncertainties observed with SHF, as uncertainties in wind speed had a much smaller impact compared to the uncertainties from MERRA-2.

7. Discussion

As seen in Figure 1 and Figure 2, CYGNSS covers nearly the entire tropical and subtropical ocean due to its orbital inclination (±35°), with a median revisit time of around 3 hours [13]. Though there are some gaps in coverage due to missing data from the Block IIF satellites [24], CYGNSS was still able to estimate the high surface heat fluxes in the northern Western Atlantic and Western Pacific Oceans (Figure 1) and those in the Southern Pacific Ocean (Figure 2). CYGNSS’s greater coverage over these parts of the world’s oceans allow it observe events such as Hurricane Florence (Figure 3). As it made landfall on the North Carolina coast, large latent heat flux values were observed by CYGNSS in its eastern half. Though Florence had weakened at this point to a Category 1, it still made a substantial impact on the Carolina coast with significant flooding as its forward speed decreased [44]. Florence’s weaker, though still hurricane strength, winds, along with a significant air-sea humidity difference, could have contributed to the higher LHF values, which in turn might have contributed to its impact on the Carolina coast.

Though CYGNSS is a tropical mission, it has the capacity to observe low-latitude extratropical cyclones that often feature large surface heat fluxes [14]. While CYGNSS only observed the equatorward side of this extratropical cyclone (Figure 4), it did observe the area where the highest fluxes are observed in ETCs. Here, the counterclockwise rotation of an ETC draws colder and drier air from the continent, which then interacts with the warm and moist ocean surface below. These air-sea differences, coupled with the high wind speeds from the cyclone, lead to significant latent and sensible heat flux values, which could affect the development and evolution of this rapidly developing ETC given the amount of energy being released from the surface and into the lower levels of this ETC. However, in both the ETC case study and Hurricane Florence, while strong surface heat flux values were observed with the CYGNSS Flux product, future studies are needed to determine how much of an impact these fluxes had on these systems.

Both the Hurricane Florence and the ETC case studies highlight some of the limits with the CYGNSS Flux product, namely with the fluxes estimated with the YSLF product. In Figure 3c,d, the LHF and SHF results using YSLF winds yield fewer samples within the storm compared to its FDS counterpart. This is due to COARE 3.5 being verified for wind speeds up to 25 m/s and any flux results above that speed limit are flagged as their quality may be suspect. Additionally, as can be seen in Figure 4c,d, the YSLF results appear noisier than the FDS fluxes. Given that time averaging is not applied to consecutive specular points for the YSLF wind speeds, it yields noisier wind speed results [16], which contributes to the noisier LHF and SHF results when using the YSLF product. Future versions of the COARE algorithm and releases of the CYGNSS wind speed product could help address some of these gaps and uncertainties in the CYGNSS Flux product.

When comparing the CYGNSS Flux product results with fluxes estimated from buoy data, one can see that while CYGNSS matches up well at lower flux values with some discrepancies, there is greater scatter at higher values, with most CYGNSS Fluxes being underestimated in comparison to buoy data at these higher values (Figure 6). While most of the fluxes between CYGNSS and buoy data fall near or along the 1-to-1 line, there is still a significant spread between the two datasets, which can be better seen in Figure 7. Here, one can see that there is a higher count of CYGNSS observations when LHF is between 100–200 W/m² and SHF below 50 W/m² (though the difference is much smaller) and a higher count of buoy observations between 0–100 W/m². At first glance, it seemed as if these differences were the result of errors and uncertainties in the CYGNSS surface wind speed product, especially at higher wind speeds where there is greater known uncertainty [40]. While an initial glance at how the differences in LHF and SHF change as wind speed increases (Figure 8) shows that CYGNSS underestimates the fluxes at higher wind speeds, a more detailed analysis indicates that MERRA-2 also contributes to these differences.

As shown in Figure 9, differences between temperature and humidity from MERRA-2 and buoy data correlated with ΔLHF and ΔSHF. While wind speeds also correlated with ΔLHF, the correlation was minimal in regard to ΔSHF. This indicates that, while there are uncertainties associated with the CYGNSS winds, the differences from MERRA-2 have a broader impact across all four of the CYGNSS Flux outputs. This can also be seen in the uncertainty analysis. While limited to just 3 days of filtered data and using the same uncertainty values across the globe for all inputs from MERRA-2, the uncertainty from these inputs on LHF and SHF were often equal or greater than the uncertainty from the CYGNSS wind speeds. This indicates that errors and uncertainties in the CYGNSS Flux product are equally coming from MERRA-2 and CYGNSS L2 wind inputs. Since one uncertainty value was used across the globe for each MERRA-2 input, future analysis could look at how uncertainties from MERRA-2 vary throughout the globe under different conditions, which in turn would affect the uncertainties in CYGNSS Flux products. Additionally, utilizing more days of data could aid in improving this simple sensitivity analysis.

8. Conclusions

By using the Level-2 surface wind speeds retrieved from CYGNSS measurements, combined with temperature and humidity estimates from MERRA-2 reanalysis data, we have been able to develop a Surface Heat Flux Product for the CYGNSS mission. The fluxes in this product were estimated using the COARE 3.5 algorithm, which has been shown in previous research to reliably estimate surface fluxes over the world’s oceans up to wind speeds of 25 m/s [18]. With proper application of quality control flags (i.e. removal of points associated with Block IIF GPS transmitters, points with RCG < 3, etc.) and the speed limit from the COARE algorithm, CYGNSS is able to produce estimates of LHF and SHF over the majority of its orbit. The CYGNSS Surface Heat Flux Product can be used to estimate the surface fluxes in weather phenomena such as tropical (Figure 3) and extratropical cyclones (Figure 4), as well as tropical convection and the general climate.

While direct flux measurements have been limited during the CYGNSS mission, we can use wind speed, temperature and humidity observed by buoys, apply those observations to the same flux algorithm to estimate LHF and SHF and compare them to the results from CYGNSS. As was seen in Figure 6, there is good agreement among all four flux products at lower flux values. Considering the entire range of values, LHF RMSD values were around 40 W/m², while SHF RMSD values were around 10 W/m² (Table 1). Overall, there was a close agreement between the CYGNSS Fluxes and the fluxes estimated from buoy data with the same algorithm, validating this product. Uncertainty did increase with increasing flux values in these comparisons, which was related to uncertainties in the CYGNSS wind speed retrievals, especially at higher wind speeds, as well as uncertainties in the inputs from MERRA-2. Future versions of the CYGNSS Surface Heat Flux Product will address these errors and will take into account future improvements of the CYGNSS L2 winds, as well as improvements in MERRA-2 or utilize a different source for temperature and humidity. Additionally, future work could look at how the CYGNSS Flux Product compares to other algorithms and surface heat flux products, similar to what has been done in previous studies [45,46].

The CYGNSS Surface Heat Flux Product described in this paper can aid the scientific community in their understanding of how heat fluxes correlate with, and possibly influence, various weather phenomena. While direct in-situ surface flux observations will continue to be the standard for ocean surface heat flux observations, the CYGNSS Surface Heat Flux Product can provide reliable coverage between observational gaps over the tropical and subtropical oceans with its high temporal and spatial frequency. The CYGNSS Flux product is now available to the public through the Physical Oceanography Distributed Active Archive Center (PO.DAAC), with data ranging from the start of the CYGNSS science mission (18 March 2017) through the present, with a time lag of about 1 month [47].

Author Contributions

J.A.C. and D.J.P. developed the CYGNSS Surface Heat Flux Product, J.A.C. was involved with the data curation, formal analysis, investigation and writing the original draft of this manuscript. D.J.P. was involved with funding acquisition and reviewing and editing the manuscript. S.A. provided the Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, Table 1 and writing in Section 5 for the product validation, as well as reviewing and editing the manuscript.

Funding

This work was supported by NASA CYGNSS Science Team Grant NNH17ZDA001N and by the CYGNSS mission under NASA Science Mission Directorate Contract NNL13AQ00C.

Acknowledgments

We would like to thank: Jim Edson at Woods Hole Oceanographic Institution and Aaron Paget of Concord University for their initial assistance with the COARE 3.5 algorithm, J. Brent Roberts at NASA Marshall for preprocessing, standardizing and providing the buoy data used in Section 5 and Catherine Naud for providing valuable feedback on the manuscript. We would also like to thank members of the CYGNSS Science Operations Center (SOC) and the Physical Oceanography Distributed Active Archive Center (PO.DAAC) for their assistance in developing the CYGNSS Surface Heat Flux Product and releasing it to the public. The comments and suggestions given by the four anonymous reviewers were helpful and greatly appreciated. This work was carried out at the Jet Propulsion Laboratory, California Institute of Technology under a contract with the National Aeronautics and Space Administration. Copyright 2019 California Institute of Technology. Government sponsorship acknowledged.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Full day observations of Latent (a,c) and Sensible (b,d) Heat Fluxes [W/m2] from Cyclone Global Navigation Satellite System (CYGNSS) (FDS version) on 1 January 2018 from 0000–1159 UTC (a,b) and 1200–2359 UTC (c,d).

Figure 2. Full day observations of Latent (a,c) and Sensible (b,d) Heat Fluxes [W/m2] from CYGNSS (FDS version) on 1 July 2018 from 0000–1159 UTC (a,b) and 1200–2359 UTC (c,d).

Figure 3. CYGNSS Surface Heat Flux Observations of Hurricane Florence on 14 September 2018. Contours of sea level pressure [black lines] at 1800 UTC with CYGNSS observations spanning ±3 hours around this time. (a) LHF with FDS winds, (b) SHF with FDS winds; (c) LHF with YSLF winds; (d) SHF with YSLF winds.

Figure 4. CYGNSS Flux Observations of an extratropical cyclone (ETC) on 4 January 2018. Contours of sea level pressure [black lines] at 1500 UTC with CYGNSS observations spanning ±3 hours from this period. (a) LHF with FDS winds, (b) SHF with FDS winds; (c) LHF with YSLF winds; (d) SHF with YSLF winds.

Figure 5. Location of buoys used for CYGNSS Flux Product validation.

Figure 6. Two dimensional-density plot of collocated CYGNSS and Buoy surface latent (a,c) and sensible (b,d) surface heat fluxes. Comparisons for both versions of CYGNSS Fluxes (top: Fully Developed Seas, bottom: Young Seas with Limited Fetch). The diagonal gray line is the 1:1 agreement.

Figure 7. Histogram of surface latent (a,c) and sensible (b,d) surface heat flux comparisons, binned at 25 W/m2, between CYGNSS and buoy data as seen in Figure 5, with buoy on the left and CYGNSS right in each plot.

Figure 8. Differences in fluxes (CYGNSS-Buoy) and how they vary with increasing wind speed observations from the buoy data. This is done for LHF (a,c) and SHF (b,d) and their respective FDS (a,b) and YSLF (c,d) products.

Figure 9. Latent (LHF) and sensible (SHF) heat flux residual dependency for each bulk component. Residual data were normalized to the range of −10 to 10. The slope of the linear fit and correlation values (star marks show statistically significant value at 95% level) are included in the bottom right of each respective plot. Flux differences (CYGNSS–buoy) for (a–c) LHF FDS, (d–f) LHF YSLF, (g–i) SHF FDS and (j–l) SHF YSLF and how they all vary with wind (ΔU). Specific humidity at 10-meters (Δqa) and surface (Δqs) is shown for LHF and temperature at the surface (ΔTa) and 10-meters (ΔTs) for SHF. The solid line represents the linear regression.

Table 1

Statistical parameters, calculated from the density scatter plots in Figure 6, between the CYGNSS and buoy flux data. RMSD: root mean square difference (CYGNSS-Buoy) [W/m²]; μ: mean bias [W/m²]; σ: standard deviation of the difference [W/m²]; N: total sample size; and r: correlation coefficient.

	Latent Heat Flux		Sensible Heat Flux
	FDS	YSLF	FDS	YSLF
RMSD	38.38	39.17	9.15	8.88
μ	3.33	3.53	−0.3	−0.37
σ	38.24	39.01	9.15	8.87
N	21,679	20,947	21,679	20,947
r	0.8	0.78	0.85	0.85

Table 2

Standard deviation of the LHF and SHF products (in W m²) when factoring in uncertainties from CYGNSS L2 winds and MERRA-2 variables. Standard deviations from the uncertainties of the individual flux components (wind, temperature, humidity) are also shown.

	Latent Heat Flux		Sensible Heat Flux
	FDS	YSLF	FDS	YSLF
Wind Speed	24.1	35.9	2.1	3.1
Surface Temperature	15.1	15.1	4.8	4.8
10m Temperature	24.2	24.2	9.7	9.6
Relative Humidity	21.3	21.2	0.2	0.2
Overall	43.9	51.4	11.2	11.4