Turbulent Heat Flux Reconstruction in the North Pacific from 1921 to 2014

Turbulent heat flux is the main passageway for air–sea interactions. However, owing to a lack of long-term observations of turbulent heat flux, it is difficult to investigate the mechanisms of coupled ocean–atmosphere variabilities, such as the Pacific Decadal Oscillation. In this study, we reconstructed the long-term turbulent heat flux in the North Pacific from 1921 to 2014 on the basis of observations in the International Comprehensive Ocean-Atmosphere Data Set–International Maritime Meteorological Archive. Sea surface temperature, air temperature, wind, and humidity were used to reconstruct the turbulent heat flux by using the Coupled Ocean– Atmosphere Response Experiment 3.5 algorithm. The modified Fisher–Tippett distribution was employed to calculate the turbulent heat flux at each grid square, and missing values were further derived on the basis of data interpolating empirical orthogonal functions. The reconstructed turbulent heat flux was shown to be in accordance with the commonly used short-term heat flux datasets. This reconstruction is further examined by comparing it with long-term data from the European Center for Medium-Range Weather Forecasts twentieth-century reanalysis (ERA-20C) and the Twentieth Century Reanalysis (20CR) dataset from the National Oceanic and Atmospheric Administration. This reconstruction displays good agreement with ERA-20C both in spatial and temporal scales but shows some differences from 20CR. These examinations show that the reconstructed turbulent heat flux can reproduce well the main features of air–sea interaction in the North Pacific, which can be used in the studies of air–sea interaction in the North Pacific on multidecadal timescales.


Introduction
Research on decadal changes in climate has attracted considerable interest during the past two decades (e.g., Kushnir 1994;Zhang et al. 1997;Mantua et al. 1997;Minobe 1997;Latif et al. 2007;Newman et al. 2016;Zhang et al. 2019). There are two main modes in the multidecadal variability of the Northern Hemisphere: the Atlantic Multidecadal Oscillation (AMO) (Delworth et al. 2007) and the Pacific Decadal Oscillation (PDO) (Mantua et al. 1997). Furthermore, many investigations have been conducted on their mechanisms (e.g., Delworth et al. 1993;Delworth and Greatbatch 2000;Newman et al. 2003;Miller and Schneider 2000). The mechanism of PDO is less understood than that of AMO (Liu 2012;Newman et al. 2016). Air-sea interactions are a fundamental problem in discussing the mechanisms of the decadalscale variability of the climate. As the main air-sea passageway, turbulent heat fluxes (latent and sensible heat fluxes) are useful in understanding the possible physical processes of air-sea interactions (e.g., Cayan 1992;Tanimoto et al. 2003;Taguchi et al. 2012;Liang and Yu 2016). However, long-term records of turbulent heat flux in observations in the North Pacific are unavailable; therefore, turbulent heat flux cannot be used to investigate air-sea interactions in multidecadal time series such as PDO. The long-term datasets of turbulent heat flux in observations are required for the analyses of the multidecadal variability of Earth's climate.
Several different datasets of turbulent heat flux are available, including the Objectively Analyzed Air-Sea Flux (OAFlux) dataset (Yu et al. 2008) from Woods Hole Oceanographic Institution, the Reanalysis Flux dataset from the National Centers for Environmental Prediction (NCEP) (Kalnay et al. 1996), the National Oceanography Centre Southampton 2.0 (NOCS2.0) Surface Flux dataset (Berry and Kent 2011), the Japanese 55-year Reanalysis (JRA-55) (Kobayashi et al. 2015;Harada et al. 2016), the Geophysical Fluid Dynamics Laboratory version 2 forcing for Common Ocean-Ice Reference Experiments (CORE2) (Large and Yeager 2009), the European Center for Medium-Range Weather Forecasts twentieth-century reanalysis (ERA-20C) (Poli et al. 2016), and the Twentieth Century Reanalysis (20CR) dataset from the National Oceanic and Atmospheric Administration (Compo et al. 2011). The OAFlux, NCEP, JRA-55, and CORE2 datasets start from the 1950s, and the NOCS2.0 dataset is available from 1973. These datasets are too short to use in multidecadal climate investigations. The ERA-20C and 20CR heat flux dataset are sufficiently long and have been used in some studies (Gao et al. 2016). Nevertheless, these two reanalysis datasets may contain biases because of model uncertainties, and some difference exists between the two. It is necessary to reconstruct a long time series of turbulent heat flux on the basis of observations. At present, the marine observations obtained by ships have long time series of variables related to turbulent heat flux. It is possible to construct a long time series for turbulent heat flux by observations from ships (Bunker 1976;Isemer and Hasse 1985;Josey et al. 1999). It is notable that observations from ships contain biases, and many works have attempted to reduce the biases (Berry et al. 2004;Berry and Kent 2011).
Currently available long-term marine observations are collated in the International Comprehensive Ocean-Atmosphere Data Set (ICOADS)-International Maritime Meteorological Archive 3.0 (Freeman et al. 2017). The ICOADS 3.0 dataset provides some airsea variables, such as sea surface temperature (SST), air temperature, wind, and humidity, thus indicating that it is possible to reconstruct the turbulent heat flux. A traditional Coupled Ocean-Atmosphere Response Experiment (COARE3.5) algorithm (Edson et al. 2013) is often used to construct the turbulent heat flux. According to this algorithm, there are two methods for calculating the turbulent heat flux. The first method is the classical approach, in which flux-related variables are gridded, and the heat flux is calculated from these gridded values (Oberhuber 1988;Berry and Kent 2009). The second method is the sampling approach, in which the heat flux is calculated at each observation site (Gulev 1997). Gulev et al. (2013) employed the latter to construct the turbulent heat flux in the North Atlantic. In the current study, we used the sampling approach to reconstruct the turbulent heat flux in the North Pacific in the domain 120°E -130°W, 25 -50°N from 1921 to 2014. Data from 2015 are also provided in ICOADS, and the version since 2015 is version 3.0.1, which has some differences from version 3.0.0; thus, data after 2014 were not considered in this research.
This paper is structured as follows: Section 2 introduces the source data and bias correction. Section 3 describes the number of observations and reconstruction of the relative humidity. Section 4 presents the reconstruction methods. Section 5 shows the testing of the quality of reconstructed data via comparisons with multiple data. Section 6 presents the examination of the reconstructed heat flux by air-sea interactions in the North Pacific. Section 7 shows the discussion and conclusions.

Source data and bias correction
The ICOADS 3.0 dataset covers the most widely observed surface marine meteorological data collected from 1662 to 2014. To reduce possible inhomogene-ities from the different types of observations, only reports from ships were selected to construct the turbulent heat flux. Several variables are involved in each marine observation. Among these variables, SST, air temperature, wind, and atmospheric humidity are required in turbulent heat flux reconstruction. However, these four variables are not always included simultaneously at each observation report in the ICOADS 3.0 dataset. For example, some observation reports include the SST and air temperature but not wind or humidity. Compared with other variables, the number of observations of atmospheric humidity is limited, particularly in the earlier decades, wherein the number of observations was too small for sampling. To supplement these data, relative humidity was first reconstructed (specific method introduced in next section), and the observations (including the reconstructed relative humidity) containing all four variables simultaneously were then selected to reconstruct the turbulent heat flux.
It has been reported that there are systematic biases in the variables obtained from marine observations (Berry and Kent 2011). Researchers have tried to reduce these errors (Josey et al. 1999), and we used the following processes to minimize the bias.
(1) For air temperature, biases exist because of the solar heating of the instruments and their surroundings. Berry et al. (2004) tried to reduce the heating biases by using a heat budget, and their method was employed in our research. Several variables, such as wind speed, cloud cover, and dew-point, are required in this method. However, there are few observations for some variables (e.g., cloud cover and dew-point) before the 1970s, thus making it difficult to estimate solar heating bias. In the current research, we employed the ERA-20CM cloud cover and dew-point reanalysis. The results based on ERA-20CM and the observations since 1970s were compared, and the comparison results show that the difference is negligible, with a root mean squares error (RMSE) less than 0.07°C (figure not shown). Thus, ERA-20CM cloud cover and dew-point reanalysis were adopted to reduce the bias before the 1970s. It was also indicated that there are biases in the night marine air temperature (Kent et al. 2013), particularly during the Second World War. Day-night difference (Kent et al. 2013) was employed to correct the night marine air temperature from 1942 to 1945. Moreover, given that the heights of the air temperature observations are different, temperature values from different heights were adjusted to a common reference height (10 m) (Kent et al. 2013).
(2) It was believed that there was a spurious upward trend in wind speed (Thomas et al. 2008), and scientists have made efforts to reduce this phenomenon. We followed the method provided by Tokinaga and Xie (2011) to remove the biases and then adjusted the measured wind speed to a common reference height (10 m). Considering that data on anemometer heights are not available before 1970, the measured wind was employed after 1970. Some anemometer height records after 1970 are missing, and a monthly mean of the 5° area gridded dataset of anemometer heights was considered the default value (Berry and Kent 2009). For the estimated wind speed, Lindau (1995) estimates were employed (Kent and Taylor 1997) to reduce bias. It was also reported that the night wind speed was underestimated because of poor sight visibility; thus, a day-night difference was added on the night wind speed to make a correction (Thomas et al. 2008;Tokinaga and Xie 2011). The estimated wind speed after 1980 was not used in our research because some reports on the estimated wind were actually measured by anemometers after 1980 (Tokinaga and Xie 2011). Therefore, the estimated wind speed employed in this paper was from 1921 to 1980.
(3) Biases also exist in SST. Different observation methods may result in different levels of bias. There are three methods for ship SST measurement. Some use a bucket to haul a sample of water to measure the SST, including the uninsulated (canvas) bucket and insulated (wooden or rubber) bucket. The water sample in the bucket may be cooled in most cases while hauling. The heat loss in a canvas bucket is larger than that in a wooden or rubber bucket. The second ship measurement is known as engine room intake (ERI) measurement. This measurement may be influenced by the depth of the observation, circumstances peculiar to each vessel, and heating of the water sample by the superstructure of the ship. The third method is measurement by using dedicated hull contact sensors. Hull sensors can provide more consistent SST measurements than ERI, but the bias still needs to be evaluated (Kent et al. 2010). We followed the method of HadSST3 (Kennedy et al. 2011) to adjust these biases.
The SST required in COARE3.5 is skin temperature. However, the SST we obtained is the bulk temperature. To transfer the bulk temperature to the skin temperature, shortwave and longwave radiation are required. COARE3.5 provides the code for transferring the bulk SST into skin temperature. The shortwave radiation and longwave radiation were set at default values in ICOADS3.5. However, these values could vary with region and time and could influence the results. We further evaluated whether this influence was sensitive. Considering the limited observations for solar radiation and thermal radiation in earlier decades, the ERA-20CM shortwave and longwave radiation were employed to calculate the heat flux. The heat flux based on default and varying values are almost the same (the RMSE of the latent and sensible heat flux is only 0.84 and 0.41 W m −2 , respectively). It is shown that the latent and sensible heat fluxes are not sensitive to the choice of radiation values.
(4) The method of Berry and Kent (2011) was employed to adjust the specific humidity. The relative humidity was first calculated by air temperature, dewpoint, and sea level pressure (SLP). Owing to the rare observations for the relative humidity in the early decades, regression was conducted to reconstruct the relative humidity. Thereafter, the specific humidity was calculated on the basis of relative humidity, air temperature, and SLP. To reduce bias from different measurements (i.e., marine screen and sling psychrometers), the specific humidity was multiplied by 0.966 for the measurement of the marine screen based on a previous work (Berry and Kent 2011). For the unknown measurement methods, 30 -70 % of the unknown methods for the measurement of specific humidity were randomly selected to be multiplied by 0.966. To further reduce the possible random errors, it was repeated for 30 times, and the ensemble mean was considered the adjusted specific humidity. After the adjustment, the residual bias uncertainty was calculated to be 0.25 g kg −1 on the basis of the differences between the observed humidity from the two different observation methods (marine screens and sling psychrometers). This result is similar to the result of Berry and Kent (2011). Figure 1a shows the mean number of reports (reports containing the four variables simultaneously) for each 5° × 5° grid square in each season. There are more than 100 reports in each grid square after 1940 ( Fig. 1a), less than 45 reports per grid square in many places, and even less than 20 reports in some areas before 1940 (Fig. 1a). Considering that there were few observations of humidity in earlier decades, the relative humidity was reconstructed using multivariable linear regression of the wind speed, air temperature, and SLP (Gulev et al. 2013) (the robustness of the relative humidity obtained through the linear regression is introduced in the next second paragraph). To obtain a better reconstruction, we tested the regression results for different time periods, and the period of 1991 -2010 was selected to identify the regression formula for the reconstruction of relative humidity.

Number of observations and reconstruction of relative humidity
After the regression, the number of observations for the reconstructed heat flux in each grid square was larger than before, particularly during 1921 -1940. The number of observations after regression was more than 20 in most locations before 1940, much larger than the original number (Fig. 1a). This is illustrated by the time series shown in Fig. 1b. For the original data (blue line), the number of observations was very small before 1950 and increased rapidly after 1950, reaching a peak in the mid-1960s, when the number of observations from ships was at a maximum. Thereafter, the number decreased, and many other types of observations became available, e.g., from buoys. By regressing the relative humidity, the number of observations increased in most decades, particularly in the 1920s and 1930s. The reason for this is that there were many reports of SST, air temperature, and wind speed, but the number of reports containing humidity is limited in those decades. After regression, the number of reports in the 1940s is the smallest compared with other decades. Regressing the relative humidity increased the number of marine observations and made it possible to calculate the turbulent heat flux.
Wind speed, air temperature, and SLP were employed to reconstruct the relative humidity. The regression model was performed for each grid and all seasons, and the result of the reconstructed relative humidity was examined. Figures 2a and 2b show the climatology of the reconstructed humidity and the difference in the climatology between the original and reconstructed data. The distribution of the humidity showed a south-north pattern, with a larger magnitude in the north and a smaller magnitude in the south (Fig. 2a); this result was the same as that for the original data (figure not shown). In most locations, the differences in climatology between the reconstructed and original relative humidity were less than ± 2.0 % (Fig. 2b), which is very small relative to the climatology. Figure 2c shows the latent heat fluxes based on the observed and reconstructed relative humidity. The RMSE is 8.76 W m −2 with a slope of 0.94, thus showing that the reconstruction of relative humidity introduced small errors. Figure 2d shows the sensible heat fluxes based on the original and reconstructed relative humidity, with an RMSE of 0.94 W m −2 . The changes of relative humidity do not affect the sensible heat flux result obviously because the relative humidity is the main factor in latent heat flux calculation but not in sensible heat flux. Figures 2e and 2f show the scatter plot of the heat flux based on the original and climatological mean relative humidity. The RMSE is 11.98 W m −2 for latent heat flux from the two methods (Fig. 2e). The RMSE is 1.22 W m −2 for sensible heat flux (Fig 2f), larger than the values shown in Figs. 2c and 2d, thus indicating that humidity reconstruction by the regression method is better than the method of using a climatological value. The left-hand column shows the original numbers, whereas the right-hand column shows the numbers after the reconstruction of the relative humidity. (b) Time series of the number of observations for each 5° × 5° grid square in different years. The blue line shows the original number, and the red line shows the number after the reconstruction of the relative humidity. It should be noted that the number in (b) is in a logarithmic scale.

Methods for the reconstruction of heat flux
4.1 COARE3.5 algorithm The COARE3.5 algorithm is widely used to calculate the turbulent heat flux. Liu et al. (1979) developed bulk aerodynamic formulas on the basis of Monin-Obukhov similarity, and the COARE algorithm is an example of this approach. There are several versions of the algorithm (Webster and Lukas 1992;Fairall et al. 1996;Brunke et al. 2002;Fairall et al. 2003), and Edson et al. (2013) updated the bulk flux algorithm to the COARE3.5 algorithm by improving the parameterization of the drag coefficients.
The formulas for the latent heat flux Q LH and the sensible heat flux Q SH are as follows: where ρ is the air mass density with the unit of [kg m −3 ], and L e the latent heat of evaporation with the unit of [J kg −1 ]. C e and C h are the turbulent exchange coefficients for the latent and sensible heat fluxes, respectively, and are functions of wind speed, detection height, and atmospheric stability. C p is the specific heat capacity of air at a constant pressure with the unit of [J kg −1 °C −1 ], and U is the scalar wind speed with the unit of [m s −1 ]. q s and q a are the specific humidity over sea water and the atmospheric specific humidity near the sea surface, respectively, with the unit of [g kg −1 ] and the unit being transferred to [kg kg −1 ] during calculations. T s is the SST, and T a air temperature with the unit of [°C]. Edson et al. (2013) reported the detailed calculations of coefficient, ρ, L e , C e , C h , and C p .

Modified Fisher-Tippett distribution
The heat flux was reconstructed by the COARE3.5 algorithm at each report and was processed to obtain values for each 5° × 5° grid square. Figure 1a shows the mean number of observations for each grid square in each season. There were large regional inhomogeneities, with fewer observations available in the central North Pacific and more observations near the coastline. The number of observations also varied greatly over time and was smaller in the early decades than in recent decades (Fig. 1b). To minimize these inhomogeneities in both space and time, we randomly selected the same number of observations in each 5° × 5° grid square, e.g., 45 (or 20) reports for each season. We calculated the heat flux for every season (DJF as winter, MAM as spring, JJA as winter, and SON as autumn) for each grid square instead of for each month as a result of the limited number of observations in the early decades. Traditionally, the mean of the 45 (20) reports would be used as the grid value, but the probability density function (PDF) of the turbulent heat flux does not follow a normal distribution (Gulev and Belyaev 2012) and cannot be used. Gulev and Belyaev (2012) studied the probability distributions of surface turbulent heat fluxes. They found that a modified Fisher-Tippett (MFT) distribution provides a good description of the statistical properties of turbulent heat fluxes. Therefore, MFT distribution was employed in this research.

a. Formulas of the MFT distribution
The MFT distribution was applied to the heat flux calculation for each 5° × 5° grid square. In each grid square, 20 (or 45) samplings were selected to calculate the grid value. The PDF of the turbulent heat flux is as follows: where α is nondimensional, and β has the dimensions of the inverse units of a variable. The two parameters can be calculated as follows: where n is the sampling number of the heat flux data (20 or 45 in each grid square), and x is the turbulent heat flux value for each sample. Given that the analytical solution cannot be obtained using Eq.
(3), we calculated the numerical solution instead. To avoid a floating overflow, the turbulent heat flux was scaled by 10 3 , and the value of β was multiplied by 10 3 . By using the values of parameters α and β, the value of the latent and sensible heat flux for each grid can be computed as follows: For each 5° × 5° grid square in each season, 20 (or 45) samplings were selected. The parameter α and β were calculated by Eq.
(3). The turbulent heat flux for each grid was then obtained using Eq. (4). By the process above, the turbulent heat flux for each 5° × 5° grid square was acquired.

b. Results of the MFT distribution
We selected 20 samples for each 5° × 5° grid square. Figure 3 plots the result of parameters α and β. Parameter α for the latent heat flux was between 2 and 3.5 in most locations in the North Pacific; this value was larger at lower latitudes and smaller at high latitudes. Parameter α for the sensible heat flux was between zero and two and was slightly larger in the western North Pacific. The magnitude and distribution off α or both latent and sensible heat fluxes are in good agreement with previously reported works (Gulev and Belyaev 2012). The magnitudes off β or both latent and sensible heat flux were between 0 and 250. The distribution off β or the latent heat flux was smaller in the southwest North Pacific and larger in the north, whereas β for the sensible heat flux was smaller in the northwest North Pacific and larger in the south. The magnitudes and distribution of β are in good agreement with previously reported work (Gulev and Belyaev 2012), thus showing that our MFT distribution is reasonable.
We also tried to randomly select 45 samples for each 5° × 5° grid square. Figures 4a and 4b show the difference in climatology between the methods with 20 and 45 samples for both the latent and sensible heat fluxes from 1950 to 2010. The difference between the two sampling methods is less than 2 W m −2 in most areas; this value is very small relative to the climatology. The largest difference is located in the low latitude and center of the North Pacific, where the numbers of observation are relatively small (Fig. 1a). A pointto-point comparison is shown as a scatter plot (Figs. 4c, d). The slope in Fig. 4c is 0.99 and has an RMSE of 5.5 W m −2 (latent heat flux), and the sensible heat flux has a slope of 1.01 and an RMSE of 2.8 W m −2 , thus showing that only a slight difference existed between the two samplings. It seems that the differences are negligible for both reconstructions on the basis of the 20 and 45 samplings. Figures 4c and 4d show the uncertainties in sampling (5.5 W m −2 for the latent heat flux and 2.8 W m −2 for the sensible heat flux).
Given that the number of reports was less than 45 in many areas in the early decades, particularly before 1950, we selected 20 reports in the entire period. There may be random errors because of the random selection of the 20 reports. To minimize these errors, we repeated the random selection 20 times, and the ensemble mean of the 20 experiments was taken as the heat flux value over the 5° × 5° grid square.

Data Interpolating Empirical Orthogonal
Functions (DINEOF) By using the processes described in the preceding sections, values were obtained for each 5° × 5° grid square. However, there are some missing data in some areas and periods. Figure 5 shows the number of missing grid squares in temporal and spatial scale. It shows that the number of missing grid squares is small after the 1950s and largest in the 1940s (Fig. 5a). The total number of grids without missing data is 146, with approximately 100 grids missing in the 1940s. The reason for the large numbers of missing data in that decade may be because many ships stopped sailing during the Second World War (Gulev et al. 2013), particularly in the North Pacific, thus leading to an absence of observations. Figure 5b shows the distribution of the amount of missing data from 1921 to 2014. The dark blue color in the east and west North Pacific indicates a small amount of missing datasets (less than 10). The missing datasets is the largest near the west seashore of the North Pacific and slightly smaller in the low latitude of the North Pacific. Considering these missing values, it is necessary to interpolate the missing heat flux. The data interpolating empirical orthogonal functions (DINEOF) (Beckers and Rixen 2003;Alvera-Azcarate et al. 2005) was used to obtain a full dataset.  DINEOF is a widely accepted method for reconstructing missing datasets (Huang et al. 2017;Hernández-Carrasco et al. 2018). It is a self-consistent, parameter-free technique that has the advantage of not needing a priori information. On the basis of EOF, the missing data can be reconstructed using the EOF modes and their amplitude time series. Compared with some other interpolation method such as optimal interpolation reconstruction, the results conducted by both methods are similar, but the computational time of DINEOF can be reduced 30 times (Alvera-Azcarate et al. 2005). In this method, the climatology is subtracted to derive the anomalies, and parts of nonmissing data are randomly removed for cross validation. Thereafter, the missing and removed data were set to zero. EOF analyses were conducted to decompose the heat flux, and the optimal number of EOF modes was identified by cross validation with randomly removed data. Finally, the missing data were replaced by EOF reconstruction. Figure 6 shows the RMSEs for different numbers of EOF modes. When the RMSE is the smallest, the corresponding EOF mode was selected to construct the heat flux. The RMSE are large for the first several modes and decreases to a relative small value as the number of modes increases. Thereafter, the RMSE grows again. For latent heat flux, the minimum RMSE occurs at the 13th EOF mode with a value of 15.51 W m −2 . The minimum RMSE for sensible heat flux also occurs at the 12th EOF mode with a value of 7.48 W m −2 . Therefore, it may bring uncertainty of 15.51 and 7.48 W m −2 for latent and sensible heat fluxes by DINEOF, respectively. On the basis of the processes above, the latent and sensible heat fluxes from 1921 to 2014 of the whole domain in the North Pacific can be acquired. To make a comparison with the reconstruction, all of these datasets were degraded to the resolution of 5° × 5°. It should be noted that a positive value indicates that heat is released from the ocean to the atmosphere in the reconstruction, OAFlux, NCEP, NOCS2.0, and JRA-55. However, the opposite is true in CORE2. For consistency, we defined the positive value when heat is released from the ocean to the atmosphere.

Comparison with heat flux datasets from recent
The climatology of the five turbulent heat flux datasets is compared in Fig. 7. For the latent heat flux, all six datasets presented a south-north pattern: the latent heat flux is smaller in the north and larger in the south, with the largest latent heat flux located over the Kuroshio extension (KE). For the sensible heat flux, an east-west pattern is shown: the sensible heat flux is smaller in the east and larger in the west, with the largest sensible heat flux located over the KE. Therefore, the reconstructed latent and sensible heat flux are consistent with the available heat flux in terms of climatology.
The correlation coefficients between the reconstructed data and the other datasets were also computed, and For the latent heat flux, the correlation coefficients between the reconstructed data and OAFlux data were > 0.4 in most areas (Fig. 8a). The largest correlation is located in the west of the North Pacific, with a value of approximately 0.8. A similar correlation pattern can be seen for the relation between the reconstructed data and NCEP, NOCS2.0, JRA-55, and CORE2 datasets (Figs. 8b -e). Figures 8a -8e show that the correlations between the reconstruction and NOCS2.0 (Fig. 8c) are the largest possibly because both the reconstructed and NOCS2.0 heat fluxes are from marine observations, whereas the others may be heat fluxes from others sources and models as an example. The correlations for sensible heat flux present similar results (Figs. 8f -j). To understand the somewhat weak correlation between the reconstruction and other heat flux data, the flux-related variables (air temperature, SST, wind speed, and specific humidity) of the reconstruction were correlated with other source data ( Table 1). The correlation between the adjusted and other source air temperature are high (larger than 0.9); these are also high for SST at approximately 0.80. The correlations are a little smaller for wind speed at approximately 0.7. The correlations between the adjusted and other source data for specific humidity vary for different data: 0.87 between the adjusted and NOCS2.0 and 0.37 between the adjusted and NCEP. It seems that uncertainties in humidity from different data are relatively large. We conclude that the reason for the somewhat weak correlation between the reconstructed and other source heat fluxes in some areas may be the result of humidity and wind. Overall, the reconstruction correlated well with other commonly used heat flux data. Figure 9 shows the time series of these turbulent heat flux data during 1974 -2006. Given that the North Pacific is large, we divided the region into four domains marked in Fig. 8e as domain A (120°E -170°W, year to year. The reconstructed data correlate well with the OAFlux, NCEP, JRA-55, NOCS2.0, and CORE2 datasets, particularly for the latent heat flux (all > 0.85). Similar results were found for the other three domains (Fig. 9). The time series of the reconstruction are relatively in good agreement with other commonly used turbulent heat flux data. The trends of the reconstruction, OAFlux, NCEP, NOCS2.0, JRA-55, and CORE2 heat flux datasets were compared during 1974 -2006 (Fig. 10). These datasets were anomalies from the seasonal average. The trends of these latent heat flux datasets show similar patterns, with positive trends in most areas of the domain, particularly in the southwest North Pacific. There are also some differences between these datasets. The largest positive trend is located in the northwest direction for the OAFlux, NCEP, and JRA-55 but is located eastward for the reconstruction and NOCS2.0. The trends of reconstruction are similar to the NOCS2.0. Both of these results are from marine observations, whereas the others are from reanalysis or a combination of observations. The sensible heat flux also shows a similar result, with the reconstructed and NOCS2.0 showing similarities. Therefore, the trends of the reconstruction are similar to other datasets, particularly NOCS2.0. The analyses above show that the reconstruction is in good agreement with other commonly used heat flux data in recent decades.

Comparison with long-term heat flux datasets
The reconstruction was further compared with longterm heat flux data  to test how the reconstruction performed on a long time period. Considering that the ERA-20C is from 1900 to 2010 and the 20CR started from 1871, the common period of 1921 -2010 was selected in performing the comparisons. The resolution for ERA-20C is 1° × 1°, and 20CR has a T62 Gaussian grid with 192 × 94 points. Both heat flux data were degraded to a resolution of 5° × 5° to perform a comparison with the reconstruction. All of these datasets were anomalies from seasonal climatology. Figure 11 shows the correlations between the recon-struction and the other two long-term heat flux data. For the latent heat flux, the correlations between the reconstruction and ERA-20C are > 0.4 in most areas (Fig. 11a); this value is larger in the south than in the north. The correlation between the reconstruction and 20CR is different: the positive correlation is located in the west and east of the North Pacific, whereas the negative correlation is located in the center of the North Pacific (Fig. 11b). The positive areas are areas with enough observations, whereas negative areas have insufficient observations (Figs. 1a,5b). In areas with sufficient observations, the reconstruction and 20CR latent heat flux show large consistency, whereas some differences existed in areas with deficient obser- vations. This is not the case in the correlation between the reconstruction and ERA-20C latent heat flux. In the latter, the correlations are high in areas regardless of the number of observations. The correlations between the reconstruction and ERA-20C sensible heat flux are larger than 0.4 in most areas; this value is a slightly smaller than that of the latent heat flux (Fig. 11c). A similar pattern can be seen in the correlations between reconstruction and 20CR sensible heat flux (Fig. 11d). The patterns in Figs. 11b and 11d are different. The correlations in the west and east North Pacific are both positive but are negative in the center of the North Pacific. According to Eq. (1), except for the parameters and wind speed, latent heat flux is determined by humidity, and sensible heat flux is determined by air temperature and SST. It is concluded that the negative correlation in the center of the North Pacific in Fig. 11c may be due to the inconsistency of the reconstructed and 20CR humidity.

Air-sea interaction in the North Pacific
The reconstructed heat flux was applied to the airsea interaction in the North Pacific to test whether it can reproduce the result of previous research. Former research shows that the atmosphere drives the ocean in the North Pacific in boreal winter in most areas (Cayan 1992;Iwasaka and Wallace 1995). When the Aleutian low has strong anomaly, enhanced wind speed and reduced air temperature and humidity will cool the ocean via turbulent heat flux (Newman et al. 2016), and Ekman transport would also amplify this pattern (Alexander and Scott 2008). However, in the subarctic frontal zone (SAFZ), specifically in the KE and Oyashio frontal zones, positive anomaly SST may release heat flux to the atmosphere (Tanimoto et al. 2003;Taguchi et al. 2012). We correlated the heat flux and SST in the North Pacific to see whether the reconstruction can reproduce the phenomenon. The SST is HadISST with a resolution of 1° × 1° from the Met Office Hadley Centre (Rayner et al. 2003). To be correlated with the turbulent heat flux, the resolution of SST was degraded to 5° × 5°. Figure 12 shows the correlation between the two variables in the North Pacific in boreal winter. The correlations are negative in most areas when the atmosphere leads the SST for a month (DJF average for turbulent [latent plus sensible] heat flux and JFM average for SST), thus indicating that the strengthening of the surface westerlies may enhance the underlying SST (Fig. 12a) (similar to Fig. 1c in Tanimoto et al. 2003). When SST leads the heat flux for a month (NDJ average for SST and DJF average for turbulent heat flux), strong positive correlations can be seen in the SAFZ (Fig. 12b). A similar result can also be seen in the lead-lag correlation between the ERA-20 turbulent heat flux and SST ( Figs. 12c, d). These patterns are similar to Fig. 2c in Tanimoto et al. (2003), thus indicating that the ocean drives the atmosphere in that area (Miller et al. 1998;Deser et al. 1999). The reproduction of the air-sea interaction in the North Pacific by the reconstructed turbulent heat flux indicates that our reconstruction is reasonable.

Conclusions and discussion
The turbulent heat flux in the North Pacific from 1921 to 2014 was reconstructed on the basis of the ICOADS 3.0 dataset by using the COARE3.5 algorithm. The SST, air temperature, wind speed, and humidity were required to calculate the heat flux data. Given that insufficient humidity data were available in the earlier decades, we first constructed the relative humidity via regression on the basis of wind speed, air temperature, and SLP. The regressed humidity was in good agreement with the actual observed humidity. The humidity was then combined with other three variables (SST, air temperature, and scalar wind) to reconstruct the turbulent heat flux according to the COARE3.5 algorithm. To minimize the inhomogeneities in the sampling, 20 reports were randomly selected over each 5° × 5° grid square for each season to calculate the grid value. The MFT distribution was used to calculate each grid value instead of the average value. To reduce the sampling errors, we repeated the experiment 20 times, and the ensemble mean was taken as the value for the grid square. However, even when using this method, there were still some missing data in some locations during some periods, particularly in the central North Pacific and in the 1940s. DINEOF was used to interpolate the missing values. An uncertainty of 17.22 W m −2 for the latent heat flux and 7.51 W m −2 for the sensible heat flux was introduced by the DINEOF reconstruction. These processes above were used to obtain the latent and sensible heat fluxes in each grid square in the North Pacific from 1921 to 2014. The reconstructed heat flux was shown to be in accordance with the commonly used heat flux datasets (OAFlux, NCEP, NOCS2.0, JRA-55, and CORE2) for recent decades. The climatology patterns of the reconstruction coincide with the other turbulent heat flux data, with the largest over the KE for both latent and sensible heat fluxes. The correlations between the reconstruction and other heat flux data are relatively high in most areas, except for some minor correlations in some areas, which may be the result of the uncertainty of the relative humidity and wind speed. The trends of the reconstruction show similar patterns to other commonly used turbulent heat flux. The reconstructed heat flux was further compared with the longterm data (ERA-20C and 20CR). The reconstructed heat flux (both latent and sensible) is highly correlated with the ERA-20C. Correlations between the reconstructed and 20CR latent heat fluxes show a pattern of positive values in the west and east of the North Pacific and negative values in the center of the North Pacific. This may be the result of the inconsistency of the reconstructed and 20CR relative humidity. The reconstructed turbulent heat flux can reproduce the airsea interaction in the North Pacific: the strengthening of the surface westerlies enhances the underlying SST in most areas, and the ocean drives the atmosphere in the SAFZ. These examinations indicate the rationality of the reconstruction. The long-term datasets of the reconstructed turbulent heat flux can be used to study the physical processes of the air-sea interaction in the North Pacific on multidecadal timescales. Future work will focus on the related physical processes, and we will try to identify which process predominates in the airsea interaction in the North Pacific, e.g., whether the atmosphere influences the ocean (Hasselmann 1976;Frankignoul and Reynolds 1983) or the ocean forces the atmosphere (Deser et al. 1999;Qiu et al. 2010). This will improve our understanding of multidecadal variability in the North Pacific.