Oceanic Evaporation: Trends and Variability

where sea surface saturated humidity is determined by sea surface temperature (SST) and salinity, and ρ is density of moist air. The transfer coefficient is dependent on the stability of the atmosphere and the sea state (Liu et al., 1979; Zeng et al., 1988). Historically, marine surface observations have provided the basis for estimating these oceanic turbulent fluxes (e.g. Bunker, 1976; Cayan, 1992; da Silva et al., 1994; Esbensen & Kushnir, 1981; Hastenrath, 1980; Hsiung, 1985; Isemer & Hasse, 1985, 1987; Josey et al., 1998; Oberhuber, 1988; Renfrew et al., 2002; Weare et al., 1981). The advent of remote sensing techniques offers means to retrieve a number of surface bulk variables. Microwave radiation interacts directly with water molecules and hence is effective in providing water vapor information. The sea surface emissivity is affected by the sea state and foam conditions, which is related to surface wind. For instance, global microwave measurement of the Special Sensor Microwave Imager (SSM/I) on board a series of Defense Meteorological Satellite Program (DMSP) satellites has been used to retrieve near-surface air humidity and winds over the ocean.


Introduction
The global water and energy cycles are strongly coupled as two essential components of earth system. They play important roles in altering the Earth's climate.
Oceanic evaporation, or sea surface latent heat flux (LHF) divided by latent heat of vaporization (L v ), is a key component of global water and energy cycle. In a bulk aerodynamic formulation, it is determined by the transfer coefficient of evaporation, C E , and bulk parameters such as surface wind speed (U), surface saturated and near-surface air specific humidity (Q s and Q a ), where sea surface saturated humidity is determined by sea surface temperature (SST) and salinity, and ρ is density of moist air. The transfer coefficient is dependent on the stability of the atmosphere and the sea state (Liu et al., 1979;Zeng et al., 1988). Historically, marine surface observations have provided the basis for estimating these oceanic turbulent fluxes (e.g. Bunker, 1976;Cayan, 1992;da Silva et al., 1994;Esbensen & Kushnir, 1981;Hastenrath, 1980;Hsiung, 1985;Isemer & Hasse, 1985Josey et al., 1998;Oberhuber, 1988;Renfrew et al., 2002;Weare et al., 1981). The advent of remote sensing techniques offers means to retrieve a number of surface bulk variables. Microwave radiation interacts directly with water molecules and hence is effective in providing water vapor information. The sea surface emissivity is affected by the sea state and foam conditions, which is related to surface wind. For instance, global microwave measurement of the Special Sensor Microwave Imager (SSM/I) on board a series of Defense Meteorological Satellite Program (DMSP) satellites has been used to retrieve near-surface air humidity and winds over the ocean.
At present there are several remote sensing products of global ocean surface latent heat flux. They include the NASA/Goddard Satellite-based Surface Turbulent Flux (GSSTF) dataset version 1 (Chou et al., 1997) and version 2 (GSSTF2, Chou et al., 2003), the Japanese Ocean Flux utilizing Remote Sensing Observations (J-OFURO) dataset (Kubota et al., 2002) and the Hamburg Ocean Atmosphere Parameters and Fluxes from Satellite (HOAPS) dataset (Grassl et al., 2000). Chiu et al. (2008) examined "trends" and variations in these global oceanic evaporation products for the period 1988-2000. They found a long-term increase in global average LHF that started around 1990 in GSSTF2. They argued that the dominant patterns may be related to an enhancement of Hadley circulation and El Niño-Southern Oscillation

Data and methodology
Earlier version of these flux products have been described elsewhere (Chiu et al., 2008). The product versions described here represent the most updated versions as of the writing of this report.

HOAPS
Detail descriptions of the latest version of HOAPS, (version 3, or HOAPS-3) are given in Andersson et al. (2010). Bulk variables are derived from SSM/I data except for the SST which is derived from the Advanced Very High Resolution Radiometer (AVHRR) Oceans Pathfinder SST product. A neural network algorithm is used to derive U. The Q a is obtained using the linear relationship of Bentamy et al. (2003). The Q s is computed from the AVHRR SST using the Magnus formula (Murray, 1967) with a constant salinity correction factor of 0.98. The near-surface air temperature (T a ) is estimated from the SST using the assumptions of 80% constant relative humidity and a constant surface-air temperature difference of 1 K. Latent and sensible heat fluxes are calculated using the Coupled Ocean-Atmosphere Response Experiment (COARE) 2.6a bulk algorithm (Fairall et al., 1996(Fairall et al., , 2003. The HOAPS-3 data sets cover the time period from July 1987 to December 2005. HOAPS-G pentad and monthly data sets with 0.5-degree resolution and HOAPS-C twice daily data set with 1-degree resolution are available at the website (http://www.hoaps.zmaw.de).

J-OFURO
The updated version of J-OFURO, (version 2, J-OFURO2) is described in Tomita et al. (2010). Bulk variables U, Q a and SST (Q s ) are determined by multi-satellite and multiple satellite sensors. U is obtained from a combination of microwave radiometers (SSM/I, AMSR-E and TMI) and scatterometers (ERS-1, ERS-2 and QuikSCAT). Q a is derived from SSM/I measurements. SST is taken from the Merged satellite and in-situ data Global Daily SST (MGDSST, Sakurai et al. 2005) analysis provided by Japanese Meteorological Agency (JMA). T a is obtained from NCEP/DOE reanalysis. COARE 3.0 bulk algorithm (Fairall et al., 2003) is used to estimate LHF and SHF. The J-OFURO2 covers the time period from January 1988 to December 2006. Daily and monthly means with 1-degree resolution are available at the website (http://dtsv.scc.u-tokai.ac.jp/j-ofuro).

GSSTF
The GSSTF2 product has daily and monthly fields with a 1°x1° resolution for July 1987-December 2000 (Chou et al., 2001), based on the method of Chou et al. (1997) with some improvements (Chou et al., 2003). The temporal and spatial resolutions of GSSTF2b are the same as those of GSSTF2, except that GSSTF2b product covers a longer period (July 1987-December 2008).
GSSTF2b dataset is processed using improved input datasets, namely the recently released NCEP SST analysis, and a uniform (across satellites) surface wind and microwave brightness temperature (TB) V6 dataset from the SSM/I produced by RSS. Table 1 summarizes characteristics of input data and parameters for HOAPS3, J-OFURO2, GSSTF2 and GSSTF2b, in that order. As we focus on LHF, only detailed descriptions and discussions on input parameters of LHF are presented.
A major improvement in the input parameters of GSSTF2b is the use of the newly released SSM/I V6 product (see discussions in http://www.ssmi.com). The SSM/I V6 product removes the spurious wind speed trends found in the Wentz/RSS SSM/I V4 wind speed retrievals. To be consistent, the SSM/I V4 total precipitable water (W) and bottom-layer precipitable water (WB) used in GSSTF2 are replaced by the corresponding SSM/I V6 products in the production of GSSTF2b. Moreover, the weekly 1° spatial resolution Optimum Interpolation (OI) SST version 1 (V1) dataset (Reynolds & Smith, 1994) used in GSSTF2 is replaced by the improved OI SST version 2 (V2) dataset. The OI SST V2 has a lower satellite bias, a new sea ice algorithm, and an improved OI analysis (Reynolds et al., 2002) resulting in a modest reduction of the satellite bias and global residual biases of roughly −0.03°C. The major improvement in the V2 analysis shows up at high latitudes where local differences between the old and new analysis can exceed 1 °C due to the application of a new sea ice algorithm. There are two GSSTF2b sets, Set1 and Set2 (Shie et al., 2010;Shie 2010). Set1 is developed using all available DMSP SSM/I sensor data. In a preliminary analysis, it was noted that there are large trends associated with LHF which are mostly attributed to the DMSP F13 and F15 satellites. Set2 was produced by excluding satellite retrievals that are judged to caused relatively large artificial trends in LHF (mostly post 2000) from Set1. Consequently Set2 is identical to Set1 before 1997 and shows a smaller trend than Set1 for the whole period, while Set1 has better spatial coverage (less missing data).  further found a drift in the Earth incidence angle (EIA) associated with the SSM/I sensors on the DMSP satellites that introduces artificial trends in the SSM/I TB data. These artificial trends introduce large changes in the boundary water (WB), which affects the Q a , and thus the LHF retrievals. An improved version, GSSTF2c, incorporating the corrected SSM/I brightness temperature, has been produced as of this writing . The retrieved WB, Q a and LHF have genuinely improved, particularly in the trends post 2000 ). An extensive study involving the GSSTF2c will be presented in a separate paper.
In this chapter, "trend" is used to indicate results from linear regression analysis and/or Empirical Mode Decomposition (EMD) for the period of study. Linear regression of the time series with time is used to detect linear trends and the significance can be estimated from the slope of the regression. EMD is based on local characteristic time scales of the data and is therefore applicable for analyzing nonlinear and non-stationary processes (Huang et al., 1998). It decomposes the time series into a finite and often small number of intrinsic mode (1996,2003) Fairall et al. functions (IMFs) of increasing time scales. The existence of a trend is dependent on the length of the dataset. If the last IMF (one with longest time scale) is monotonically increasing or decreasing, a trend is indicated. The EMD is a more stringent test for significance of "trends." Non-seasonal variability is examined using Empirical Orthogonal Function (EOF) analysis. Non-seasonal data are obtained by subtracting the monthly climatology of the study period from the monthly data. EOF analysis decomposes a spatio-temporal dataset into a series of orthogonal spatial EOF patterns and the associated time series (also called principal components). The test proposed by North et al. (1982) is used to judge the EOFs to see if they are significant and distinct. To examine the significance of each EOF, the logarithm of the variance explained by the EOF is plotted against its EOF number. The variance explained by the nth EOF is given by where λ n is the nth eigenvalue and N is the number of time samples. Linear regression between the logarithms of λ n vs n is computed. EOFs above the regression line are judged to be significant. Figure 1 shows the area weighted global (60°N-60°S) average LHF of HOAPS3, J-OFURO2, GSSTF2, and Set1 and Set2 of GSSTF2b. The merged OAFLUX product is included for comparison. Visual inspection of the data products does not show large missing gaps in the spatial distribution of the products. All satellite products-HOAPS3, J-OFURO2 and all GSSTF datasets show increases while there is no obvious trend in OAFLUX. All products have similar global means in the early period 1988-1991. The means of GSSTF2b (both Set1 and Set2) are generally lower than that of GSSTF2 which is the highest among all products. All products show a dip in 1991 and an increase afterwards. The dip is clearly evident in HOAPS3. GSSTF2b Set1 and Set2 are identical up to 1998 after which they diverge, but tend to come close again after 2006. The zonal annual averages of LHF are depicted in Figure 2. The general features of the zonal means are quite similar among these products: they showed maxima in the subtropics, minima at the poles and relative minima at the equator. The subtropical maxima in the southern hemisphere are slighter higher than those in the northern hemisphere for the same product. While GSSTF2 is the highest among these estimates, the GSSTF2b Set2 is slightly lower than HOAPS3 but higher than J-OFURO2 and OAFLUX at their maxima. Poleward of 30°, the GSSTF zonal means are generally higher than the other products.
To map out the geographic differences, Figure

EMD analysis of global average LHF
We performed EMD analyses for both GSSTF2 and GSSTF2b for the periods 1988-2000 and 1988−2008, respectively. Figure 4 shows the last IMFs from EMD analysis of GSSTF2b global average (60°N-60°S) LHF for the period 1988-2008, along with those of HOAPS3, J-OFURO2 and OAFLUX. EMD analysis of GSSTF2 global average LHF for the period 1988-2000 (Chiu et al., 2008;Xing, 2006) showed an increase starting around 1990 before the Mt Pinatubo eruption event in 1991. The last IMF of GSSTF2b Set2 and HOAPS3 show a dip around 1990 while J-OFURO2 and GSSTF2b Set1 show monotonic increases. From EMD analysis, a "trend" cannot be ascertained for the whole period for GSSTF2b Set2 and HOAPS3.

GSSTF2b Set1 and Set2
In this subsection, we analyze the difference between GSSTF2b Set1 and Set2. The linear trends of GSSTF2b LHF for the entire period 1988-2008 are depicted in Figure 5. The magnitudes of the trends are reduced when compared to GSSTF2 (see Chiu et al., 2008) and Set1 shows larger trends than Set2, as anticipated.
Trends of the zonal averages are shown in Figure 6. Both Set1 and Set2 show similar trend patterns while the magnitude in Set2 is much reduced. Maximum increasing trends are located in the subtropics, with maxima at 30°S and 35°N, poleward of the latitude of maximum zonal evaporation, while the trend at the thermal equator (~5°N) shows a minimum.  . Unit in W m −2 decade −1 . Contours give the trends above 95% confidence level.  Figure 7 shows the spatial distribution of trends in U, DQ and Q a in Set 1 and Set2. It can be seen that the trend pattern in U is almost identical in both Set1 and Set2. Major difference is found in the DQ field, the magnitudes of the trends in Set1 are larger than that in Set2. In the equatorial eastern Pacific, in portion of the South China Sea, and the Bay of Bengal, Set2 actually shows a decreasing but non-significant trend. For the Q a trend, Set1 shows large areas of decreases. While the patterns are similar for Set2, the magnitudes are much reduced, and in some areas, such as the eastern North Pacific and the bulk of the North Atlantic, the decreasing trends actually reverse to increasing trends.
LHF is a product of the surface wind (U) and the humidity difference (DQ=Q s -Q a ) if we assume the variation in C E is small (Equation 1). Judging from the change pattern of U and DQ they are essentially decoupled. Equation (1) can be integrated globally to get where δx represent the change in the quantity x, δx/x represent fractional changes in x, and the over-bar x represents global average of x. For GSSTF2 (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000), the terms in equation (2) is approximately 17%, 6%, and 11%, in that order (Xing, 2006). Most of the increase in DQ was attributed to increase in s Q and decrease in a Q . Figure 8 shows the time series of global average U, DQ and Q a for GSSTF2b Set1 and Set2, respectively. It clearly indicates the divergence of Set1 and Set2 in late 1997. The large difference in LHF between Set1 and Set2 is mostly attributed to DQ, which is due to the higher Q a in Set2. The difference in U between the datasets is small. There is a large decrease in Q a at the end of 2008. Units are in m s −1 decade −1 for U and in g kg -1 decade -1 for DQ and Q a . Contours give the trends above 95% confidence level. Table 3 shows the changes in LHF in Set1 and Set2 of GSSTF2b (1988GSSTF2b ( −2008 and the associated changes in U, DQ and the changes in Q s and Q a . The change in DQ contribute most to the change in LHF for both Set1 and Set2, while the changes in DQ is due both to an increase in Q s and decrease in Q a . The difference in the change of LHF between Set1 (23.1%) and Set2 (15.5%) is mostly attributed to DQ (20.0% vs. 12.3%) and changes in Q a (-0.51 g kg -1 vs. -0.24 g kg -1 ). It is clear that the impact of DMSP F13 is the introduction of a much lower Q a , thus affecting DQ and ultimately LHF.

EOF analyses and teleconnections
Empirical Orthogonal Function (EOF) analyses are performed on LHF of both GSSTF2b Set1 and Set2 for the period 1988-2008. Monthly means for the entire period are first removed from the data to form the non-seasonal dataset. The first, second and third EOF of Set1 Fig. 9. Spatial patterns of the first (top) and second (bottom) EOF of GSSTF2b LHF for the period 1988-2008. The variances explained are 8.6% (11.5%) and 4.3% (4.3%) for EOF1 and EOF2 for Set2 (Set1), respectively. explains 11.5%, 4.3%, and 3.4% of the total variance, and the corresponding variance explained for Set2 are 8.6%, 4.3% and 3.4%. Figure 9 shows the spatial patterns of the first (EOF1) and second (EOF2) EOF for Set1 and Set2, respectively. Their associated time series, accompanied by a rescaled Southern Oscillation Index (SOI), are presented in Figure 10. The general patterns of EOF1 of Set1 and Set2 are very similar, with large weights in subtropical Indian Ocean, the dry zone in the eastern tropical South Pacific and South Atlantic. They also bear striking resemblance to the first EOF pattern computed from GSSTF2 which shows weights of the same sign everywhere (Chiu et al., 2008). However, in GSSTF2b, there are regions of opposite sign in the South China Sea area in both Set1 and Set2, with higher weights in Set2 in the centers of highs. The EOF2 patterns for Set1 and Set2 are almost identical. Similarity is also noted for the EOF2 pattern derived from GSSTF2 (Chiu et al., 2008).
The time series of EOF1 also show a slight increasing trend, suggesting an increase of LHF over most of the global ocean. The associated time series of EOF2 have a significant correlation of 0.73 (0.72 for Set1) with SOI, reaffirming that EOF2 -ENSO events association.
The third EOF pattern (not shown) is characterized by a negative -positive -negative (-+ -) zonal changes centering around 10−20°N, 20−30°N and north of 40°N in the north Pacific and a (+ -+) centering at 20°N, 30°N and 50°N in the North Atlantic. The pattern correlations for EOF1, EOF2 and EOF3 between Set1 and Set2 are 0.08, 0.49, and 0.39, and the corresponding temporal correlation for the associated time series are 0.92, 0.97, and 0.83, in that order.
EOF3 is reminiscent of the North Atlantic Oscillation (NAO) and the North Pacific Oscillation (NPO) patterns (Walker & Bliss, 1932;Wallace & Gultzer, 1981, Hurrell et al., 2003. These teleconnection patterns are further discussed in terms of an atmospheric annular mode, or Arctic Oscillation (AO), showing the opposition between subtropical highs and the polar lows (Wallace, 2000;Deser, 2000;Aubaum et al., 2001). We compute the correlations between an AO index with the time series of EOF3 and found no significant correlation with Set2, however a correlation of 0.32 is found for EOF3 Set1, significant at 95% level.

Summary and discussion
Four satellite based sea surface latent heat flux (LHF) products, HOSAPS3, J-OFURO2, GSSTF2, and GSSTF2b (Set1 and Set2) and a merged analysis OAFLUX are compared. Linear trend analysis of all satellite based products show large increasing trend, with GSSTF2 the largest, followed by GSSTF2b Set1, J-OFURO2, HOAPS3, and GSSTF2b Set2. OAFLUX exhibits the lowest linear increasing trend. Most of the satellite products used SSM/I as input. Small drifts in the SSM/I brightness temperature (TB) associated with changes in Earth incidence angle (EIA) was noted in most of the SSM/I data . Because of the sensitivity of the boundary layer water (WB) to the TB, these small drifts can introduce artificial trends in bulk quantities such as Q a . A second data set, GSSTF2b Set2, which excludes satellite retrievals that were judged to introduce these biases, was introduced. The most-excluded satellite data are SSM/I onboard the DMSP F13 and F15 satellites. The new set, GSSTF2b Set2, was found to have a much reduced increasing trend, the magnitude of which is comparable to HOAPS3 and J-OFURO2 for the period of overlap. To account for the drift in the EIA, a new version of GSSTF, GSSTF2c, that takes account of the correction in EIA, has been completed as of this writing, and will be officially released to the public via NASA/GES DISC by the end of October 2011 . These trend issues will be revisited after its release.
Empirical Mode Decomposition (EMD) analyses, which are designed for examining nonstationary non-homogeneous time series, are performed on the global LHF. The last IMF of GSSTF2b Set1 shows a monotonic increase indicating the existence of a trend in this period.
The corresponding IMF of the other data products do not show monotonic increases, hence trends cannot be ascertained.
To examine the attribution of the increase in GSSTF2b Set1 and Set2, the linear trends in both the surface wind (U) and surface humidity difference (DQ) are computed. There is no significant difference between the trend patterns in the wind field. However, large difference in the DQ trend is noted. The DQ trend difference is attributed to a reduction in the negative trend in surface air humidity (Q a ) in Set2.
A major difference between GSSTF2 and GSSTF2b is the use of RSS SSM/I V4 for GSSTF2 and RSS SSM/I V6 for GSSTF2b. The changes in LHF, U, and DQ are 16%, 6%, and 11% for GSSTF2 (Xing, 2006). The corresponding changes are 23%, 3%, and 20% for GSSTF2b Set1 and 16%, 3%, and 12% for Set2 (Table 3). The use of RSS SSM/I V6 products reduces the wind trend from 6% to 3% but increases the DQ trend from 11% to 20% for Set1. The exclusion of F13 and F15 data reduces the LHF trend to 16%, mostly due to a reduction in the DQ trend to 12% (from 20%) with U changes remain at 3%.
Interannual variability is examined using EOF analyses. The first three significant nonseasonal EOF patterns are similar, and they explaining 10.5%, 4.3% and 3.4% for Set1 and 8.6%, 4.3% and 3.4% for Set2, respectively. The first EOF pattern of GSSTF2 for 1998-2000, with opposite changes between the equatorial eastern Pacific and the subtropics in the Pacific and Indian ocean, may be indicative of an enhance Hadley circulation (Chiu & Xing, 2004;Chiu et al., 2008). Observations also indicate large decadal variability in the Hadley Circulation (Wielicki et al., 2002;Cess & Udelhofen, 2003;Chen et al., 2002;Mitas and Clements, 2005).
This seesaw pattern is much reduced in the EOF1 pattern in both GSSTF2b Set1 and Set2 of 1998-2008, which may indicate a reduction, change of phase, or mixing of the signal with the trend in GSSTF2b. The contribution to the total variance is smaller for Set2, which excluded DMSP datasets that contains large long-term trends introduced by drifts in the Earth incidence angle in the SSM/I sensors. The difference in the fraction of variance explained in Set1 and Set2 is attributed to the artificial trend in F13 and F15 and the EIA drift effect.
Examination of the trends of the zonal means show that the latitude of maximum increase, situated in the subtropics, is found poleward of the LHF maximum in the tropic. This pattern is consistent with the expansion of the Hadley Circulation associated with global warming as predicted in climate models (Lu et al., 2007).
The EOF2 patterns of Set1 and Set2 are almost identical, both contributed to 4.3% of the variance of the dataset. The association with the El Nino/Southern Oscillation phenomena is corroborated by a high correlation between their time series and an index of the Southern Oscillation (SOI). The patterns for EOF3 and their associated time series are also similar, indicating that both GSSTF2b Set1 and Set2 are useful for examining interannual variability.

Acknowledgment
This study is supported by the MEaSUREs Program of NASA Science Mission Directorate-Earth Science Division. The authors are especially grateful to their program manager M. Maiden and program scientist J. Entin for their valuable supports of this research.