Estimation of Surface Air Specific Humidity and Air – Sea Latent Heat Flux Using FY-3 C Microwave Observations

Latent heat flux (LHF) plays an important role in the global hydrological cycle and is therefore necessary to understand global climate variability. It has been reported that the near-surface specific humidity is a major source of error for satellite-derived LHF. Here, a new empirical model relating multichannel brightness temperatures (TB) obtained from the Fengyun-3 (FY-3C) microwave radiometer and sea surface air specific humidity (Qa) is proposed. It is based on the relationship between TB, Qa, sea surface temperature (SST), and water vapor scale height. Compared with in situ data, the new satellite-derived Qa and LHF both exhibit better statistical results than previous estimates. For Qa, the bias, root mean square difference (RMSD), and the correlation coefficient (R2) between satellite and buoy in the mid-latitude region are 0.08 g/kg, 1.76 g/kg, and 0.92, respectively. For LHF, the bias, RMSD, and R2 are 2.40 W/m2, 34.24 W/m2, and 0.87, respectively. The satellite-derived Qa are also compared with National Oceanic and Atmospheric Administration (NOAA) Cooperative Institute for Research in Environmental Sciences (CIRES) humidity datasets, with a bias, RMSD, and R2 of 0.02 g/kg, 1.02 g/kg, and 0.98, respectively. The proposed method can also be extended in the future to observations from other space-borne microwave radiometers.


Introduction
Air-sea latent heat flux (LHF), combined with momentum flux and sensible heat flux, represents an important aspect of the atmosphere-ocean interaction, energy budget, water cycle variability, and global climate systems [1].In addition, inland lakes influence greatly on surface energy budgets and have shown considerable changes worldwide [2,3].This is particularly true in high-latitude regions, where anthropogenic climate change is predicted to be exceptionally rapid [4].LHF can be estimated by the bulk aerodynamic formula that includes sea surface temperature (SST), surface winds, air temperature (T a ), and near-surface air specific humidity (Q a ).These flux-related variables can be obtained from three major sources: In situ measurements, satellite observations, and numerical weather prediction models.Several previous studies have used in situ measurements from ships or buoys to estimate LHF [5,6].Since the temporal and spatial resolution of ship measurements is limited, large uncertainties are still present in the estimation of LHF.The atmospheric reanalysis products can be strongly influenced by variations in the type and amount of data assimilated [7].
Depending on the data fusion techniques, blended global surface heat flux datasets have been derived from multi-satellite information.These mainly include the Hamburg Ocean-Atmosphere Parameters and Fluxes from Satellite Data (HOAPS-3; [8,9]), Japanese Ocean Flux Data Sets with Use of Remote Sensing Observations (J-OFURO; [10]), Goddard Satellite-Based Surface Turbulent Flux (GSSTF-3; [11]), and French Research Institute for Exploitation of the Sea (IFREMER) turbulent flux products [12].Compared with in situ and reanalysis flux products, satellite observations provide much higher spatial resolution and a homogeneous time series of atmospheric state variables for surface flux computation.However, a large portion of the errors in satellite-derived LHF is found to be associated with uncertainties in the estimation of Q a [13].
Several authors have investigated the estimation of Q a from satellite measurements.Using the column water vapor (W) as an independent parameter is a common way to determine the Q a from satellite measurements [14,15].This Q a −W relationship is widely used in LHF-retrieval algorithms.However, errors in Q a originating from the relationship might result in large errors in estimating LHF, especially on hourly and daily scales.Schulz et al. (1993) developed a method to estimate the bottom-layer (500 m height) precipitable water from the brightness temperature [16].Thereafter, several similar studies based on Schulz's model have been proposed.For example, Schlüssel et al. (1995) (hereafter SC95) directly regressed the remotely sensed brightness temperature and Q a to avoid error propagation [17].Later, this approach was updated using improved training data [12,18].The difference in the frequencies used for Q a retrieval in each algorithm is considered to be a cause of statistical difference between each product.High-frequency channels, 85-89 GHz, had not been used for Q a retrieval in previous studies because the polarized radiation was considered to have no significant information about the boundary layer [16].Iwasaki et al. (2012) (hereafter IW12) concluded that the use of the brightness temperature determined by the 85 GHz polarized radiation can considerably reduce the Q a bias, particularly in the wet regions [19].Information on the vertical profile of humidity is also useful for improving the estimation accuracy of Q a [20,21].Moreover, methods using artificial intelligence such as neural networks and genetic algorithm have also been proposed [22][23][24].More recently, Tomita et al. (2018) (hereafter TO18) developed a method to estimate the Q a using vertical water profile information and W [25].However, the TO18 algorithm ignored the effect of SST, which is highly correlated with W, on Q a .
In this study, we developed an improved method to estimate Q a based on the TO18 algorithm.Vertical water profiles and other meteorological parameters (SST and W) were used to improve retrieval accuracy.The remaining parts of this paper are organized as follows: Datasets used in this study are introduced in Section 2. Analysis procedures and retrieval updates of Q a are described in Section 3. The results and validation are presented in Section 4. The discussion and conclusions are included in Sections 5 and 6, respectively.

Satellite Observations
The brightness temperature (T B ) was obtained by the Microwave Radiation Imager (MWRI) onboard the new generation polar-orbiting meteorological satellite of China (Fengyun-3C, FY-3C), launched in September 2013.The MWRI instrument is a microwave conical-scanning imager following on the heritage of the Special Sensor Microwave Imager/Sounder (SSMI/S), Advanced Microwave Scanning Radiometer for EOS (AMSR-E), and the AMSR-2 instruments.MWRI has vertically and horizontally polarized radiation for 10.65, 18.7, 23.8, 36.5, and 89 GHz [26] (hereafter referred as 10 V/H, 19 V/H, 23 V/H, 37 V/H, and 89 V/H respectively; the suffixes V and H mean vertical and horizontal polarization, respectively).In this study, data from ten common channels were used; the correlation between multiple channels is shown in Table A1.After data preprocessing, the instantaneous T B data from January to December 2014 were used for algorithm development.The main characteristics of the MWRI are shown in Table 1.

In Situ Measurements
Two kinds of in situ measurements were used as training data to derive the Q a regression formula.One was the International Comprehensive Ocean-Atmosphere Data Set Release 3.0 (ICOADS 3.0) [27].The other was from mooring buoy networks: The National Data Buoy Center (NDBC) buoys off the U.S. Atlantic, Pacific, and Gulf coasts maintained by the National Oceanic and Atmospheric Administration (NOAA); the European offshore data acquisition system (ODAS) buoys in the eastern Atlantic, maintained by the UK Met Office and Météo-France (MFUK); and the Tropical Atmosphere Ocean (TAO) buoys located in the tropical Pacific Ocean, maintained by NOAA's Pacific Marine Environment Laboratory (PMEL).
Data from January to December 2014 that included sea surface temperature, relative humidity (dew point temperature), air temperature, and pressure were used in this study.The Q a was computed with the saturation vapor by using a pressure formula [28].We removed the outlier values of Q a , which were above the upper quartile plus 1.5 × inter-quartile range (IQR) or below the lower quartile minus 1.5 × IQR [29].After carrying out quality control for Q a , the Coupled Ocean-Atmosphere Response Experiment version 3.0 (COARE3.0)model was used to adjust the Q a data to standard height (10 m) using the metadata on the type of ship and the height of the meteorological sensors [30].To match up in situ data with the satellite observations, we set criteria that the temporal and spatial differences had to be less than 30 min and 25 km, respectively.Finally, in total 23,200 pairs of ship/satellite collocated data were divided randomly into two equal groups (Samples 1 and 2).Data in Sample 1 were used for the development of the new algorithm, while Sample 2 was used for an additional evaluation of the retrieval algorithms as independent data.

NOAA CIRES
The instantaneous FY-3C Q a values are compared with NOAA CIRES Multi-Satellite Humidity.This product combines Special Sensor Microwave Imager/Special Sensor Microwave Imager Sounder (SSMI/SSMIS) and Advanced Microwave Sounding Unit (AMSU) satellite observations to create a new satellite humidity dataset.The retrievals were developed using research vessel data from the Shipboard Automated Meteorological and Oceanographic System (SAMOS) initiative and validated using the ICOADS dataset [31].The NOAA data are available through the Center for Ocean-Atmospheric Prediction Studies (COAPS) Marine Data Center (https://mdc.coaps.fsu.edu/data).

ERA-Interim Reanalysis Data
The ERA-Interim is a global atmospheric reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) [32].ERA-Interim products cover the period from 1979, and are continuously updated in real time.In this study, humidity, SST, and W were obtained from the ERA-Interim daily and gridded products, which have been used for theoretical research and validation.This data can be download from the website: http://apps.ecmwf.int/datasets/.The time series is from January 2014 to December 2014.

Methodology
There are different satellite remote sensing algorithms for air-sea specific humidity and LHF, but these algorithms do not fully consider the effects of SST, W, and vertical moisture structure on Q a .In this study, a methodological framework was proposed to estimate Q a and LHF from FY-3C observations and meteorological auxiliary data (see Figure 1).In this study, a methodological framework was proposed to estimate  and LHF from FY-3C observations and meteorological auxiliary data (see Figure 1).

Relationship between SST and W
Kanemaru and Masunaga (2013) illustrated that the water vapor scale height ( ) is a good indicator of the vertical moisture gradient between the boundary layer and the free troposphere, where W is the column water vapor (kg/m 2 ),  is the density of dry air (1.2 kg/m 3 ), and  is the surface water vapor mixing ratio (g/kg) [33].Here, values of W and  were derived from ERA-Interim [32].Tomita et al. (2018) presented the value of  influencing the relationship between  and the brightness temperature.They proposed a method using additional data for  as a proxy for the vertical moisture gradient to retrieve  beyond previous linear multi-regressions.The influence of W on  is considered in their method.
Kanemaru and Masunaga (2013) also found that the relationship between SST and W changes depending on the values of SST [33].W shows contrasting sensitivities to SST depending on the SST.The SST-W relation resembles the C-C (Clausius-Clapeyron) relation [34,35] when SST is lower than 27 °C.For higher SST, the SST-W relation does not follow the Clausius-Clapeyron relation because the influence of W on  is larger than the SST contribution.Hence, the SST should be considered when we study the relationship between W and  .

Relationship between SST and W
Kanemaru and Masunaga (2013) illustrated that the water vapor scale height (H v ) is a good indicator of the vertical moisture gradient between the boundary layer and the free troposphere, where W is the column water vapor (kg/m 2 ), ρ a is the density of dry air (1.2 kg/m 3 ), and q v is the surface water vapor mixing ratio (g/kg) [33].Here, values of W and q v were derived from ERA-Interim [32].Tomita et al. (2018) presented the value of H v influencing the relationship between Q a and the brightness temperature.They proposed a method using additional data for H v as a proxy for the vertical moisture gradient to retrieve Q a beyond previous linear multi-regressions.The influence of W on Q a is considered in their method.
Remote Sens. 2019, 11, 466 5 of 19 Kanemaru and Masunaga (2013) also found that the relationship between SST and W changes depending on the values of SST [33].W shows contrasting sensitivities to SST depending on the SST.The SST-W relation resembles the C-C (Clausius-Clapeyron) relation [34,35] when SST is lower than 27 • C. For higher SST, the SST-W relation does not follow the Clausius-Clapeyron relation because the influence of W on H v is larger than the SST contribution.Hence, the SST should be considered when we study the relationship between W and Q a .
Figure 2a-e shows the relationship between the obtained SST and W with the change in H v values.The color depth represents the value of q v .Five groups of H v (≤1300 m, 1300-1800 m, 1800-2300 m, 2300-2800 m, and >2800 m) were temporally averaged for the year 2014, using the data from 60 • S to 60 • N. Lower H v corresponded to relatively lower W values.For example, when H v ≤ 1300 m, the maximum value of W was less than 35 kg/m 2 , while the maximum value reached 85 kg/m 2 if H v > 2800 m.This may be because lower H v corresponds to atmospheric columns containing moisture in relatively lower layers, and higher H v corresponds to columns containing moisture up to relatively higher layers [25].
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 20 Figure 2a-e shows the relationship between the obtained SST and W with the change in  values.The color depth represents the value of  .Five groups of  (≤1300 m, 1300-1800 m, 1800-2300 m, 2300-2800 m, and >2800 m) were temporally averaged for the year 2014, using the data from 60° S to 60° N. Lower  corresponded to relatively lower W values.For example, when  ≤ 1300 m, the maximum value of W was less than 35 kg/m 2 , while the maximum value reached 85 kg/m 2 if  > 2800 m.This may be because lower  corresponds to atmospheric columns containing moisture in relatively lower layers, and higher  corresponds to columns containing moisture up to relatively higher layers [25].Furthermore, compared with higher  values, the increase ratio of / was smaller at lower  values (see Figure 2).The ratio changed with SST, the thick line shows a cubic regression curve calculated from SST and W. Figure 2f shows the correlation between the increase ratio and the SST at each  range.When the SST was lower than 27 °C, the increase ratio had less change with the increase of SST at lower  ranges.When the SST was near 27 °C, the ratio of / was significantly increased (especially at higher  ranges), and then the increase ratio was reduced.The influence of the SST-W relation on  will be further discussed in Section 3.2.

Channel Sensitivity
Figures 3 and 4 show the relationship between  and the brightness temperature Furthermore, compared with higher H v values, the increase ratio of dW/dSST was smaller at lower H v values (see Figure 2).The ratio changed with SST, the thick line shows a cubic regression curve calculated from SST and W. Figure 2f shows the correlation between the increase ratio and the SST at each H v range.When the SST was lower than 27 • C, the increase ratio had less change with the increase of SST at lower H v ranges.When the SST was near 27 • C, the ratio of dW/dSST was significantly increased (especially at higher H v ranges), and then the increase ratio was reduced.The influence of the SST-W relation on Q a will be further discussed in Section 3.2.

Channel Sensitivity
Figures 3 and 4 show the relationship between Q a and the brightness temperature (T B,23V and T B,89H ) at five H v ranges.The brightness temperature dataset is from the FY-3C MWRI.The color depth indicates the density of data and the black line was calculated using all H v data.There was a basic positive linear relationship between brightness temperature and Q a , and the H v values influenced this relationship.When the H v was small (Figure 3a-c, under 2300 m; Figure 4a,b, under 1800 m), the relationship appeared linear.However, the relationship between Q a and brightness temperature appeared to curve at higher H v ranges (Figure 3d and e  The basic relationship between  and the brightness temperature can be understood from the semi-linear relationship between W and the brightness temperature. , and  , are well related to columnar water vapor [19,36].When the value of W is large, the atmospheric column contains more water vapor [15].This results in a higher brightness temperature due to higher microwave radiation from the increased amount of water vapor, so the basic relationship appears linear. To further explore the effect of SST on the above phenomena, we divided the data shown in Figure 3 into four groups by SST. Figure 5 shows the relationship between brightness temperature (T , ), Q , and SST at each  range.Red lines were calculated for each SST interval, and black lines were calculated using all data of each  interval.Figure 5 demonstrates that the distribution of points is influenced by SST.In the SST ≤ 10 °C range, the slope of the scatter distribution was lower than the slope of all scatter data, and their difference was the largest in this interval.As the The basic relationship between Q a and the brightness temperature can be understood from the semi-linear relationship between W and the brightness temperature.T B,23 and T B,89 are well related to columnar water vapor [19,36].When the value of W is large, the atmospheric column contains more water vapor [15].This results in a higher brightness temperature due to higher microwave radiation from the increased amount of water vapor, so the basic relationship appears linear.
To further explore the effect of SST on the above phenomena, we divided the data shown in Figure 3 into four groups by SST. Figure 5 shows the relationship between brightness temperature (T B,23V ), Q a , and SST at each H v range.Red lines were calculated for each SST interval, and black lines were calculated using all data of each H v interval.Figure 5 demonstrates that the distribution of points is influenced by SST.In the SST ≤ 10 • C range, the slope of the scatter distribution was lower than the slope of all scatter data, and their difference was the largest in this interval.As the temperature increased, the difference between them decreased.Except for H v ≤ 1300 m, the range of 20-27 • C showed minimal difference, and then the difference increased again.This phenomenon may be due to two reasons: (1) Difference in the amount of data between different temperature ranges has an impact on the results.For example, the ratio of the number of scatters in the range of 20-27 °C to the total number of scatters was the largest, which had a greater contribution to the slope of the fitting.In contrast, for the range of 0-10 °C and >27 °C the ratio was small and corresponded to a small contribution to the overall slope.(2) The linear relationship between  and  was different at each SST range, which may be affected by the SST-W relation mentioned above.The difference in the effect of the SST increment on W further caused a change in the relationship between the  and  .
In the high frequency (89 H GHz) channel, we did the same analysis as shown in Figure 6.Compared with 23 V GHz, the scatters exhibited a more curvilinear relationship in the higher  ranges.There was a significant difference in the relationship between  and  in different SST ranges.We further used the observation sensitivity (  / ) to analyze this phenomenon quantitatively, and the results are shown in Figure 7.This phenomenon may be due to two reasons: (1) Difference in the amount of data between different temperature ranges has an impact on the results.For example, the ratio of the number of scatters in the range of 20-27 • C to the total number of scatters was the largest, which had a greater contribution to the slope of the fitting.In contrast, for the range of 0-10 • C and >27 • C the ratio was small and corresponded to a small contribution to the overall slope.(2) The linear relationship between Q a and T B was different at each SST range, which may be affected by the SST-W relation mentioned above.The difference in the effect of the SST increment on W further caused a change in the relationship between the Q a and T B .
In the high frequency (89 H GHz) channel, we did the same analysis as shown in Figure 6.Compared with 23 V GHz, the scatters exhibited a more curvilinear relationship in the higher H v ranges.There was a significant difference in the relationship between Q a and T B in different SST ranges.
We further used the observation sensitivity (dT B /dQ a ) to analyze this phenomenon quantitatively, and the results are shown in Figure 7.   Figure 7 shows that the observation sensitivity varies with SST on the same  range, and the relationship is affected by  values.For the 23 V GHz channel, when  ≤ 2300 m, the maximum value of observation sensitivity appeared in the mid-SST range, and there was minimal difference between different SST ranges.When  > 2300 m, observation sensitivity gradually decreased with the increase of SST, which resulted in a curve distribution of scatters.In the 89 H GHz channel, there was a similar trend.When  ≤ 1300 m, the maximum value of observation sensitivity appeared at 10-20 °C.Observation sensitivity decreased as SST increased when  was larger than 1300 m.In addition, Figure 7 shows that the observation sensitivity value of 23 V GHz is lower than 89 H GHz channel as a whole.Furthermore, the variation amplitude between different SST ranges was lower Each subfigure shows a linear regression black line calculated using all data of the same H v range; the red line was calculated using the corresponding SST range data.
Figure 7 shows that the observation sensitivity varies with SST on the same H v range, and the relationship is affected by H v values.For the 23 V GHz channel, when H v ≤ 2300 m, the maximum value of observation sensitivity appeared in the mid-SST range, and there was minimal difference between different SST ranges.When H v > 2300 m, observation sensitivity gradually decreased with the increase of SST, which resulted in a curve distribution of scatters.In the 89 H GHz channel, there was a similar trend.When H v ≤ 1300 m, the maximum value of observation sensitivity appeared at 10-20 • C. Observation sensitivity decreased as SST increased when H v was larger than 1300 m.In addition, Figure 7 shows that the observation sensitivity value of 23 V GHz is lower than 89 H GHz channel as a whole.Furthermore, the variation amplitude between different SST ranges was lower than 89 H GHz. This difference has an effect on the relationship between Q a and T B .Since the observation sensitivity is an inverse value of dQ a /dT B , the smaller variation of observation sensitivity in different SST ranges generates a closer relationship between Q a and T B .Hence, the relationship appears linear at 23 V GHz radiation but appears curved at 89 H GHz. Higher H v values exhibit scatter points of both channels in the form of a curve.A similar result was found for the 23 H GHz and 89 V GHz observation channels, but was unclear for 10, 19, and 37 GHz channels.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 20 than 89 H GHz. This difference has an effect on the relationship between  and  .Since the observation sensitivity is an inverse value of  / , the smaller variation of observation sensitivity in different SST ranges generates a closer relationship between  and  .Hence, the relationship appears linear at 23 V GHz radiation but appears curved at 89 H GHz. Higher  values exhibit scatter points of both channels in the form of a curve.A similar result was found for the 23 H GHz and 89 V GHz observation channels, but was unclear for 10, 19, and 37 GHz channels.In general, it was found that observation sensitivity within the same  range changed depending on the value of SST.This trend may contribute to the curved relationship between  and  , particularly at higher  values (Figures 5 and 6).This result indicates that SST is an important parameter for  estimates and cannot be ignored.Previous algorithms based only on direct regression between  and  might produce estimation errors for some ranges of SST.

Developed Algorithm
SST has been considered as an important near-surface meteorological parameter in some artificial intelligence algorithms when retrieving  and LHF [23,37].From the relationship discussed in Sections 3.1 and 3.2, we found clear connections between in situ  and  , and W was also useful because of its relationship with  .In this study, an improved model was proposed based on the TO18 method.We suggested using all the brightness temperature data in order to obtain a more accurate estimation of the  (for more details of channel correlation and its impact on the regression please see Section 5 and Table A1).Moreover, because the relationship between  and  appeared curved at higher  ranges, the quadratic terms for  , / and  , / were used for inversion.Ultimately, the new regression equation had the following form: where  denotes the specific air humidity at 10 m in g/kg,  is the brightness temperature in Kelvin, SST is in °C,  are regression coefficients, and the p-values for these retrieved coefficients are listed in Table 2.
The p-values for regression coefficients indicate the significance in estimating the  .A small pvalue (<0.05) indicates that the regression is significant.Overall, we found that all input parameters were significant in the regression model (see Table 2).But there were some exceptions, which were related to different  intervals.For example, when  ≤ 1300 m, the p-values of  and  were larger than 0.05.In comparison, when 1300 m <  ≤ 3300 m, all p-values were less than 0.05, In general, it was found that observation sensitivity within the same H v range changed depending on the value of SST.This trend may contribute to the curved relationship between Q a and T B , particularly at higher H v values (Figures 5 and 6).This result indicates that SST is an important parameter for Q a estimates and cannot be ignored.Previous algorithms based only on direct regression between T B and Q a might produce estimation errors for some ranges of SST.

Developed Algorithm
SST has been considered as an important near-surface meteorological parameter in some artificial intelligence algorithms when retrieving Q a and LHF [23,37].From the relationship discussed in Sections 3.1 and 3.2, we found clear connections between in situ Q a and T B , and W was also useful because of its relationship with Q a .In this study, an improved model was proposed based on the TO18 method.We suggested using all the brightness temperature data in order to obtain a more accurate estimation of the Q a (for more details of channel correlation and its impact on the regression please see Section 5 and Table A1).Moreover, because the relationship between Q a and T B appeared curved at higher H v ranges, the quadratic terms for T B,23V/H and T B,89V/H were used for inversion.Ultimately, the new regression equation had the following form: where Q a denotes the specific air humidity at 10 m in g/kg, T B is the brightness temperature in Kelvin, SST is in • C, c 0-15 are regression coefficients, and the p-values for these retrieved coefficients are listed in Table 2.
The p-values for regression coefficients indicate the significance in estimating the Q a .A small p-value (<0.05) indicates that the regression is significant.Overall, we found that all input parameters were significant in the regression model (see Table 2).But there were some exceptions, which were related to different H v intervals.For example, when H v ≤ 1300 m, the p-values of c 3 and c 14 were larger than 0.05.In comparison, when 1300 m < H v ≤ 3300 m, all p-values were less than 0.05, representing very significant parameters in improved models.For c 14 , the quadratic coefficient of T B,89H in extreme H v values (≤1300 m or >3300 m) was larger than 0.05.These demonstrate that in extremely large or very small H v , the term was weakly correlated with retrieved results.In this study, we deleted the terms with p-values greater than 0.05 when obtaining the retrieved coefficients.Finally, the values of coefficients are listed in Table A3.The LHF is widely estimated using the COARE3.0bulk aerodynamic formula [30], which is suitable for both satellite and in situ data.The LHF can be expressed as follows: where C e is the turbulent exchange coefficient for latent heat flux, L e is the latent heat of evaporation, and Q s is surface humidity.U and Q a are the wind speed and specific air humidity relative to the sea surface at the height of 10 m.In this study, Q a was estimated using satellite data and U was derived from FY-3C MWRI orbit products.Q s was calculated from SST assuming saturation at the surface, where a multiplier factor of 0.98 was used when considering the reduction in vapor pressure caused by a typical salinity of 34 psu.C e was computed by the COARE3.0algorithm, which was mainly influenced by conditions of light wind, stable stratification, sea spray, and other sea-state information.In addition, FY-3C MWRI orbit products of rain rate were used as input data to estimate LHF.Detailed information of these products can be obtained from the website: http://www.nsmc.org.cn/en/NSMC/Home/Index.html.

Air Specific Humidity
In this section, satellite-derived Q a fields were compared with in situ observed Q a .We also compared the proposed model with three other existing empirical algorithms (SC95, IW12, and TO18).Details of these four Q a retrieval algorithms are given in Table 3. Figure 8 shows scatter plots of the comparison of in situ observed Q a with four satellite-derived Q a values (the Sample 2 data).The proposed model has good performance, with a root mean square difference (RMSD) of 1.53 g/kg, and a correlation coefficient (R 2 ) of 0.94.TO18 has an RMSD and R 2 of the order of 1.95 g/kg and 0.92, respectively.SC95 has an RMSD and R 2 of 2.88 g/kg, and 0.86, respectively.IW12 performed better (RMSD, 2.44 g/kg; R 2 , 0.89) than SC95 due to the use of the high-frequency 89 GHz channel.Figure 8 shows scatter plots of the comparison of in situ observed  with four satellitederived  values (the Sample 2 data).The proposed model has good performance, with a root mean square difference (RMSD) of 1.53 g/kg, and a correlation coefficient (R 2 ) of 0.94.TO18 has an RMSD and R 2 of the order of 1.95 g/kg and 0.92, respectively.SC95 has an RMSD and R 2 of 2.88 g/kg, and 0.86, respectively.IW12 performed better (RMSD, 2.44 g/kg; R 2 , 0.89) than SC95 due to the use of the high-frequency 89 GHz channel.Figure 9 shows the zonal averages of the bias (satellite derived Q a minus in situ Q a ), RMSD, and standard deviation difference (SDD) for the Sample 2 data.For the proposed model, the RMSD are less than 2.0 g/kg and biases are less than 1 g/kg over most of the latitudes.The bias (Figure 9a) for all algorithms tends to show positive values, especially in low and mid-latitude regions.This means that all four algorithms tend to overestimate Q a .Compared with the other three algorithms, the proposed model's Q a shows the smallest bias.The accuracy of satellite-derived Q a was significantly improved by using SST. Figure 9 shows the zonal averages of the bias (satellite derived  minus in situ  ), RMSD, and standard deviation difference (SDD) for the Sample 2 data.For the proposed model, the RMSD are less than 2.0 g/kg and biases are less than 1 g/kg over most of the latitudes.The bias (Figure 9a) for all algorithms tends to show positive values, especially in low and mid-latitude regions.This means that all four algorithms tend to overestimate  .Compared with the other three algorithms, the proposed model's  shows the smallest bias.The accuracy of satellite-derived  was significantly improved by using SST.

Latent Heat Flux
To assess the improvement resulting from the proposed  model on satellite-derived LHF, the satellite-derived LHF was compared with in situ LHF, which was calculated from in situ measurements (i.e., Sample 2 data).Table 4 provides some statistical results of these comparisons.The statistics were averaged in the high (45°-60°N, 45°-60°S), mid (15°-45°N, 15°-45°S), and low (15°S-15°N) latitudinal regions.
Compared with the other three existing algorithms, the proposed method had good performance in estimating LHF.There were significant statistical improvements of LHF due to these improvements of  .For example, in mid-latitude regions, the bias, RMSD, and R 2 for SC95 method were 0.66 g/kg, 3.16 g/kg, and 0.81, respectively, while those of the proposed model were 0.08 g/kg, 1.76 g/kg and 0.92, respectively.Corresponding parameters of LHF (SC95) were 14.10 W/m 2 , 51.96

Latent Heat Flux
To assess the improvement resulting from the proposed Q a model on satellite-derived LHF, the satellite-derived LHF was compared with in situ LHF, which was calculated from in situ measurements (i.e., Sample 2 data).Table 4 provides some statistical results of these comparisons.The statistics were averaged in the high (45  Compared with the other three existing algorithms, the proposed method had good performance in estimating LHF.There were significant statistical improvements of LHF due to these improvements of Q a .For example, in mid-latitude regions, the bias, RMSD, and R 2 for SC95 method were 0.66 g/kg, 3.16 g/kg, and 0.81, respectively, while those of the proposed model were 0.08 g/kg, 1.76 g/kg and 0.92, respectively.Corresponding parameters of LHF (SC95) were 14.10 W/m 2 , 51.96 W/m 2 , and 0.51, while those of the proposed model were 2.40 W/m 2 , 34.24 W/m 2 , and 0.87, respectively.The differences between satellite and in situ data were dependent on zonal locations, which were mainly influenced by the accuracy of Q a .The satellite-derived Q a fields were also compared with NOAA CIRES Multi-Satellite Humidity.To match up two Q a datasets, we set criteria that the temporal and spatial differences had to be less than 180 min and 25 km, respectively.Figure 10 shows the two datasets covering the common area on 6 October 2014.
W/m 2 , and 0.51, while those of the proposed model were 2.40 W/m 2 , 34.24 W/m 2 , and 0.87, respectively.The differences between satellite and in situ data were dependent on zonal locations, which were mainly influenced by the accuracy of  .

𝑄 Comparison with NOAA CIRES Datasets
The satellite-derived  fields were also compared with NOAA CIRES Multi-Satellite Humidity.To match up two  datasets, we set criteria that the temporal and spatial differences had to be less than 180 min and 25 km, respectively.Figure 10 shows the two datasets covering the common area on 6 October, 2014.    Figure 11 shows the results of the comparison between NOAA CIRES Multi-Satellite Humidity and FY-3C humidity (total 693,836 collocated data) on 6 October, 2014.Compared with the Multi-Satellite  dataset, the  derived from FY-3C satellite had a bias, RMSD, and R 2 of 0.02 g/kg, 1.02 g/kg, and 0.98, respectively.

Discussion
In this study, we investigated the observation sensitivity ( / ) variation with the SST at different  ranges (as described in Section 3.2).This phenomenon indicates that previous  algorithms based only on direct regression between brightness temperature and  might not be accurate enough for some SST and humidity conditions.
Hence, we introduced the SST as a key input parameter into  estimation (Equation ( 2)).Radiations from 10 microwave channels of the FY-3C satellite were all used in the proposed algorithm.It is well known that  with the same frequency but different polarizations, has a strong correlation.The Pearson correlation coefficients also clearly show that there are strong dependencies between some channels used in Equation (2) (Table A1).To better show the importance of each indicator in the proposed model, a stepwise regression was used (Table A2).The results of the stepwise regression show that when the number of parameters used in the regression model was larger than six, the improvements of R 2 and residual mean square of the regression model were very small.However, it was also shown that including more channels could improve the proposed linear model, even if the newly added channel was highly correlated to existing channels.Besides, we regressed the model coefficients according to different  ranges.In order to obtain the best fitness of the regression model, we chose to include all channels in the proposed model.Moreover, we decided to use quadratic terms for  , / and  , / because the relationship between  and brightness temperature appears curved at higher  ranges.These measures ensure a superior output of our model for different humidity conditions in the global oceans.

Discussion
In this study, we investigated the observation sensitivity (dT B /dQ a ) variation with the SST at different H v ranges (as described in Section 3.2).This phenomenon indicates that previous Q a algorithms based only on direct regression between brightness temperature and Q a might not be accurate enough for some SST and humidity conditions.
Hence, we introduced the SST as a key input parameter into Q a estimation (Equation ( 2)).Radiations from 10 microwave channels of the FY-3C satellite were all used in the proposed algorithm.It is well known that T B with the same frequency but different polarizations, has a strong correlation.The Pearson correlation coefficients also clearly show that there are strong dependencies between some channels used in Equation (2) (Table A1).To better show the importance of each indicator in the proposed model, a stepwise regression was used (Table A2).The results of the stepwise regression show that when the number of parameters used in the regression model was larger than six, the improvements of R 2 and residual mean square of the regression model were very small.However, it was also shown that including more channels could improve the proposed linear model, even if the newly added channel was highly correlated to existing channels.Besides, we regressed the model coefficients according to different H v ranges.In order to obtain the best fitness of the regression model, we chose to include all channels in the proposed model.Moreover, we decided to use quadratic terms for T B,23V/H and T B,89V/H because the relationship between Q a and brightness temperature appears curved at higher H v ranges.These measures ensure a superior output of our model for different humidity conditions in the global oceans.
It is worth mentioning that the original SC95, IW12, and TO18 algorithms were derived from different training data and may not be well suited for the FY-3C MWRI instrument.Thus, we derived new regression coefficients for three algorithms before using them to compare with the proposed method.Figure 8 shows that the proposed method has the lowest RMSD (1.53 g/kg) and the highest R 2 (0.94).In addition, satellite-derived Q a were also compared with the NOAA CIRES Multi-Satellite humidity dataset, their bias, RMSD, and R 2 were 0.02 g/kg, 1.02 g/kg, and 0.98, respectively.These comparison results indicate that adding SST as an input parameter can improve Q a and LHF estimation.There are regional differences of the outputs of the proposed model.The RMSD of Q a in mid-latitude regions (15 • -45 • N, 15 • -45 • S) was 1.76 g/kg, and the RMSD of high (45 • -60 • N, 45 • -60 • S) and low (15 • S-15 • N) latitude regions were 1.11 and 1.51 g/kg, respectively.This difference may result from the SST-W relation.Kanemaru and Masunaga (2013) concluded that the W is closely related with SST and dependent on regional position [33].
We obtained SST data from the FY-3C MWRI orbit products and H v data from atmospheric reanalysis.Therefore, any potential error and dependency for either variable may affect the accuracy of Q a .Furthermore, the training and testing dataset of in situ measurements was collected from the global sea area between 60 • S and 60 • N.This means that the proposed Q a model is also only valid for these regions.Although we believe that the form of this algorithm is universal in calculating global specific humidity, the retrieved coefficients are only available for FY-3C MWRI.For other space-borne radiometers, these coefficients need to be recalculated.

Conclusions
Latent heat flux represents an important aspect of the atmosphere-ocean interaction.Additionally, satellite-derived specific humidity has been a major source of error for the estimation of LHF.In order to improve the estimation accuracy of specific humidity, we proposed a new Q a model based on FY-3C microwave brightness temperature observations.We found the relationship of observation sensitivity (dT B /dQ a ) changed depending on the range of SST, especially at larger H v , so SST was added into the estimation model.This improved model gives better estimates of Q a and LHF when compared with three existing methods (SC95, IW12, and TO18).The satellite-derived Q a were also compared with in situ data, and their bias, RMSD, and R 2 were 0.17 g/kg, 1.53 g/kg, and 0.94, respectively.The satellite Q a were also compared with NOAA CIRES Q a datasets, with a bias, RMSD, and R 2 of 0.02 g/kg, 1.02 g/kg, and 0.98, respectively.The differences between satellite LHF and in situ LHF were dependent on zonal locations: In high latitude regions (45 • -60 • N, 45 • -60 • S), the bias, RMSD, and R 2 were 12.89 W/m 2 , 23.44 W/m 2 , and 0.72, respectively; in mid-latitude regions (15 • -45 • N, 15 • -45 • S), the bias, RMSD, and R 2 were 2.4 W/m 2 , 34.24 W/m 2 , and 0.87, respectively; in low latitude regions (15 • S-15 • N), the bias, RMSD, and R 2 were 4.87 W/m 2 , 29.05 W/m 2 , and 0.82, respectively, which were mainly influenced by the accuracy of Q a (see Table 4).
As shown in the study, SST as an additional input parameter greatly improved the estimation accuracy of Q a and LHF for the FY-3C satellite observation.The innovative aspect of the algorithm was based on the finding of a relationship between SST, W, and water vapor profile information, which provided a theoretical basis for satellite-derived Q a and LHF.Furthermore, SST and W were directly used as input parameters for estimating Q a .Both variables may introduce unknown redundant information to model outputs.Future work should explore the essence of the SST-W relationship and improve the simulation method of the vertical moisture structure.

Figure 1 .
Figure 1.Methodological framework to estimate global air-sea  and latent heat flux (LHF) using FY-3C satellite observations.

Figure 1 .
Figure 1.Methodological framework to estimate global air-sea Q a and latent heat flux (LHF) using FY-3C satellite observations.

Figure 2 .
Figure 2. (a)-(e) The relationship between sea surface temperature (SST), column water vapor (W), and surface water vapor mixing ratio  (g/kg, color) on different ranges of water vapor scale height  (m).Each figure shows a cubic regression curve calculated from SST and W. (f) SST compared to increase ratio (W/SST, kg/m 2 •°C) over five  ranges, calculated from daily data and averaged from the range of 60° S to 60° N.

Figure 2 .
Figure 2. (a-e) The relationship between sea surface temperature (SST), column water vapor (W), and surface water vapor mixing ratio q v (g/kg, color) on different ranges of water vapor scale height H v (m).Each figure shows a cubic regression curve calculated from SST and W. (f) SST compared to increase ratio (W/SST, kg/m 2 • • C) over five H v ranges, calculated from daily data and averaged from the range of 60 • S to 60 • N.
Figures3 and 4show the relationship between Q a and the brightness temperature (T B,23V and T B,89H ) at five H v ranges.The brightness temperature dataset is from the FY-3C MWRI.The color depth indicates the density of data and the black line was calculated using all H v data.There was a basic positive linear relationship between brightness temperature and Q a , and the H v values influenced this relationship.When the H v was small (Figure3a-c, under 2300 m; Figure4a,b, under 1800 m), the relationship appeared linear.However, the relationship between Q a and brightness temperature appeared to curve at higher H v ranges (Figure3d and e, over 2300 m; Figure 4c-e, over 1800 m), which is similar to the results of Tomita et al. (2018) [25].Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 20 influenced this relationship.When the  was small (Figure 3a-c, under 2300 m; Figure 4a and b, under 1800 m), the relationship appeared linear.However, the relationship between  and brightness temperature appeared to curve at higher  ranges (Figure 3d and e, over 2300 m; Figure 4c-e, over 1800 m), which is similar to the results of Tomita et al. (2018) [25].

Figure 3 .
Figure 3. Scatter diagrams showing the relationship between  , ,  , and  , using the Sample 1 data.(a) For  ≤ 1300 m; (b) for 1300 m <  ≤ 1800 m; (c) for 1800 m <  ≤ 2300 m; (d) for 2300 m <  ≤ 2800 m; (e) for  > 2800 m; and (f) for all  ranges.Each figure shows a linear regression line calculated using all  range data.Color depth indicates the density of data.

Figure 3 .
Figure 3. Scatter diagrams showing the relationship between T B,23V , Q a , and H v , using the Sample 1 data.(a) For H v ≤ 1300 m; (b) for 1300 m < H v ≤ 1800 m; (c) for 1800 m < H v ≤ 2300 m; (d) for 2300 m < H v ≤ 2800 m; (e) for H v > 2800 m; and (f) for all H v ranges.Each figure shows a linear regression line calculated using all H v range data.Color depth indicates the density of data.

Figure 4 .
Figure 4. Scatter diagrams showing the relationship between  , ,  , and  , using the Sample 1 data.(a) For  ≤ 1300 m; (b) for 1300 m <  ≤ 1800 m; (c) for 1800 m <  ≤ 2300 m; (d) for 2300 m <  ≤ 2800 m; (e) for  > 2800 m; and (f) for all  ranges.Each figure shows a linear regression line calculated using all  range data.Color depth indicates the density of data.

Figure 4 .
Figure 4. Scatter diagrams showing the relationship between T B,89H , Q a , and H v , using the Sample 1 data.(a) For H v ≤ 1300 m; (b) for 1300 m < H v ≤ 1800 m; (c) for 1800 m < H v ≤ 2300 m; (d) for 2300 m < H v ≤ 2800 m; (e) for H v > 2800 m; and (f) for all H v ranges.Each figure shows a linear regression line calculated using all H v range data.Color depth indicates the density of data.

Figure 5 .
Figure 5.The relationship between T B,23V , Q a , and SST.(a-e) For H v ≤ 1300 m; (f-j) for 1300 m < H v ≤ 1800 m; (k-o) for 1800 m < H v ≤ 2300 m; (p-t) for 2300 m < H v ≤ 2800 m; and (u-y) for H v > 2800 m.Each subfigure shows a linear regression black line calculated using all data of the same H v range; the red line was calculated using the corresponding SST range data.

Figure 6 .
Figure 6.The relationship between T B,89H , Q a , and SST.(a-e) For H v ≤ 1300 m; (f-j) for 1300 m < H v ≤ 1800 m; (k-o) for 1800 m < H v ≤ 2300 m; (p-t) for 2300 m < H v ≤ 2800 m; and (u-y) for H v > 2800 m.Each subfigure shows a linear regression black line calculated using all data of the same H v range; the red line was calculated using the corresponding SST range data.

Figure 7 .
Figure 7. SST compared with observation sensitivity from a linear fitting for each SST range.(a) For 23 V GHz and (b) for 89 H GHz.

Figure 7 .
Figure 7. SST compared with observation sensitivity from a linear fitting for each SST range.(a) For 23 V GHz and (b) for 89 H GHz.

Figure 8 .
Figure 8. Comparisons between in situ observed and four satellite-derived  values (g/kg).(a) For SC95; (b) for IW12; (c) for TO18; and (d) for the proposed model.

Figure 8 .
Figure 8. Comparisons between in situ observed and four satellite-derived Q a values (g/kg).(a) For SC95; (b) for IW12; (c) for TO18; and (d) for the proposed model.

20 Figure 9 .
Figure 9. Zonal distributions of the comparative statistics between in situ observed and satellitederived  values for four algorithms, including (a) averaged bias, (b) root mean square difference (RMSD), and standard deviation difference (SDD, plotted as red lines).Statistics were averaged for each 2° bias.

Figure 9 .
Figure 9. Zonal distributions of the comparative statistics between in situ observed and satellite-derived Q a values for four algorithms, including (a) averaged bias, (b) root mean square difference (RMSD), and standard deviation difference (SDD, plotted as red lines).Statistics were averaged for each 2 • bias.

Figure 10 .
Figure 10.The common coverage area of NOAA CIRES Multi-Satellite Humidity (dark blue strip) and FY-3C humidity (colored strip) on 6 October, 2014.Color depth indicates the value of specific air humidity.(a) Ascending orbit and (b) descending orbit.

Figure 10 .
Figure 10.The common coverage area of NOAA CIRES Multi-Satellite Humidity (dark blue strip) and FY-3C humidity (colored strip) on 6 October 2014.Color depth indicates the value of specific air humidity.(a) Ascending orbit and (b) descending orbit.

Figure 11
Figure11shows the results of the comparison between NOAA CIRES Multi-Satellite Humidity and FY-3C humidity (total 693,836 collocated data) on 6 October 2014.Compared with the Multi-Satellite Q a dataset, the Q a derived from FY-3C satellite had a bias, RMSD, and R 2 of 0.02 g/kg, 1.02 g/kg, and 0.98, respectively.

Figure 11 .
Figure 11.Results of the comparison between NOAA CIRES Multi-Satellite Humidity and FY-3C humidity on 6 October, 2014.Color depth indicates the density of data.

Figure 11 .
Figure 11.Results of the comparison between NOAA CIRES Multi-Satellite Humidity and FY-3C humidity on 6 October 2014.Color depth indicates the density of data.
* V: Vertical polarization and H: Horizontal polarization.

Table 2 .
p-Values of retrieved coefficients for different water vapor scale heights H v (m).

Table 3 .
Channels and meteorological parameters for each algorithm.

Table 3 .
Channels and meteorological parameters for each algorithm.

Table 4 .
Comparisons between satellite and in situ data for the four algorithms.
a Comparison with NOAA CIRES Datasets

Table 4 .
Comparisons between satellite and in situ data for the four algorithms.