A New Empirical Estimation Scheme for Daily Net Radiation at the Ocean Surface

: Ocean surface net radiation (R n ) is signiﬁcant in research on the Earth’s heat balance systems, air–sea interactions, and other applications. However, there have been few studies on R n until now. Based on radiative and meteorological measurements collected from 66 globally distributed moored buoys, it was found that R n was dominated by downward shortwave radiation (R ↓ g ) when the length ratio of daytime (LRD) was greater than 0.4 but dominated by downward longwave radiation (R ↓ l ) for the other cases (LRD ≤ 0.4). Therefore, an empirical scheme that includes two conditional models named Case 1 (LRD > 0.4) utilizing R ↓ g as a major input and Case 2 (LRD ≤ 0.4) utilizing R ↓ l as a major input for R n estimation was successfully developed. After validation against in situ R n , the performance of the empirical scheme was satisfactory with an overall R 2 value of 0.972, an RMSE of 9.768 Wm − 2 , and a bias of − 0.092 Wm − 2 . Speciﬁcally, the accuracies of the two conditional models were also very good, with RMSEs of 9.805 and 2.824 Wm − 2 and biases of − 0.095 and 0.346 Wm − 2 for the Case 1 and Case 2 models, respectively. However, due to the limited number of available samples, the performances of these new models were poor in coastal and high-latitude areas, and the models did not work when the LRD was too small (i.e., LRD < 0.3). Overall, the newly developed empirical scheme for R n estimation has strong potential to be widely used in practical use because of its simple format and high accuracy.


Introduction
Ocean surface all-wave net radiation (R n ) is the sum of the downward and upward shortwave and longwave radiation at the ocean surface. Mathematically [1]: where R ng is the shortwave net radiation (Wm −2 ), R nl is the longwave net radiation (Wm −2 ), R ↓ g is the downward shortwave radiation (Wm −2 ), R ↑ g is the upward shortwave radiation (Wm −2 ), R ↓ l is the downward longwave radiation (Wm −2 ), and R ↑ l is the upward longwave radiation (Wm −2 ). Downward is defined as positive in this study.
The ocean occupies 71% of the Earth's surface area, and it is the largest heat-storage body in the Earth's climate system because seawater has a large specific heat capacity to absorb and store most of the solar insolation. Ocean heat flux, which is composed of the ocean surface, radiative heat flux (R n ) and turbulent heat flux (i.e., latent heat flux and sensible heat flux) [2], plays a vital role in controlling the stability of the Earth's climate system by regulating the heat balance of the Earth system through frequently exchanging heat flux with the atmosphere [3][4][5][6][7][8]. In addition, ocean heat flux was found to be closely related to the Earth's energy imbalance (EEI) [9]. The EEI is usually reflected as a larger incident solar radiation compared with the outgoing longwave radiation of the Earth's system [3,10]; therefore, it is usually used to characterize the climate state of the Earth and global warming caused by human activities [11]. Hence, it is of great significance to obtain an accurate ocean heat flux estimates. However, as an important component of the ocean heat flux, the ocean surface R n has received less attention than turbulent fluxes, although it is also critical as the input for most physical models for long-term weather and climate prediction [12]. Therefore, accurate ocean surface R n is very meaningful for various applications.
In contrast to the land surface, R n is not routinely measured for the ocean surface; instead, the downward radiative components (R ↓ g and R ↓ l ) are frequently measured but only at very few buoy sites. Because of the paucity of measurements, the ocean surface R n is usually taken out of simple climatology or parameterizations [13]. The alternative way to obtain the ocean surface R n is from various types of products, including the satellite-based products generated from radiative transfer models (i.e., Clouds and the Earth's Radiant Energy System Synoptic Radiation Fluxes and Clouds <CERES-SYN> [14] and International Satellite Cloud Climatology Project Flux Dataset <ISCCP-FD> [15]), model reanalysis products generated from weather numerical models (i.e., ERA-Interim [16]), ship-based products generated from parameterization formulas (i.e., NOCS Flux Dataset v2.0<NOCSv2.0> [17]), and reconstructed products by combining different products with interpolation methods (i.e., a third-generation Japanese ocean-flux dataset using remote sensing observations <J-OFURO3> [2] and Objectively Analyzed Air-sea Fluxes<OAFlux> [18]). However, several studies have pointed out that large discrepancies exist among these products in ocean surface R n or radiative components, and these differences cannot be ignored [13]. Hence, it is desirable to develop a relatively simple method to estimate the ocean surface R n with high accuracy for practical applications.
The traditional way to calculate R n is to add its components (shortwave and longwave radiation) estimated through physical or empirical models separately, but it is very difficult to estimate all components accurately; therefore, the accuracy of the final obtained R n would be very likely influenced by error propagation and accumulation. According to previous studies [19,20], R n is mainly determined by shortwave radiation for most cases; hence, many empirical models were successfully developed to calculate R n directly from shortwave radiation and other ancillary information, especially on the land surface [1], and even the latest released Global Land Surface Satellite (GLASS) R n product, whose accuracy was validated to be comparable to or even better than most existing products, was generated by using this kind of algorithm [21]. Similarly, Polavarapu analyzed the relationship between the ocean surface R n and R ↓ g using survey observations collected from the Global Atmospheric Research Program-Atlantic Tropical Experiment (GATE) ship and found a close linear relationship between the two components [22], based on which a linear regression model was established to estimate the ocean surface R n . Although the validation results were satisfactory on the GATE site, the robustness and global applicability of this model is questionable because of the limited samples used, and it is known that the performance of an empirical model usually relies heavily on the quality and comprehensive representation of the observations. Fortunately, an increasing number of radiative measurement platforms have been built on the global ocean surface in recent decades, increasing awareness of the significance of the ocean surface radiative flux [23,24]; therefore, there have already been several studies on ocean surface radiative component estimation, such as ocean surface downward shortwave [25][26][27] and longwave [28][29][30] measurements, but very few studies have focused on ocean surface R n . Therefore, this is the right time to extend Polavarapu's study for ocean surface R n estimation with these increased buoy measurements.
The objective of this study is to develop a new scheme for estimating the daily ocean surface R n empirically in a simple way based on measurements collected from more than 60 global moored buoy sites. Then, the performance of this new scheme is evaluated and analyzed carefully. The organization of this paper is as follows. The data and their preprocessing are introduced in Section 2. Section 3 presents the development of the new scheme for ocean surface R n estimation. The validation results and analysis of the model are provided in Section 4. The conclusions are presented in Section 5.

Data and Preprocessing
The primary data used for modeling and validation in this study are the radiative and meteorological measurements collected from 66 globally distributed, moored buoy sites from five networks/projects. In addition, other ancillary information was also used for modeling, including the ocean surface broadband albedo, broadband emissivity, clearness index (CI), and daytime duration (Dd). Table 1 summarizes all the variables used in this study and their sources. In summary, there were 53,460 daily samples in total. Specifically, for each site, 70% of the samples were randomly selected for R n estimation model development (37,378 in total), while the other 30% (16,082 in total) were used for independent validation. Furthermore, one remotely sensed product, Clouds and the Earth's Radiant Energy System Synoptic Radiation Fluxes and Clouds (CERES SYN1deg_Ed4A), was also used for comparison in this study. More details are given below.

Measurements from Moored Buoy Sites
Radiative and meteorological observations from 66 moored buoy sites from 1988 to 2018 were collected and used in this study. Figure 1 shows the spatial distribution of the 66 moored buoy sites in five networks/projects, most of which were distributed in low-and mid-latitude oceans. However, only very few moored buoys were sparsely distributed in the mid-high latitude oceans in the Southern Hemisphere, which might be attributed to the presence of fewer "hot spots" for air-sea interactions in this region [31]. Fortunately, more buoys have been considered to be deployed at seas in the Southern Hemisphere, such as at the Brazil Current (55 • -41 • W, 48 • -28 • S) [32], the Eastern Australian Current (150 • -165 • E,  Table 2 lists the detailed information about the five observation networks to which the moored buoy sites belong. The Upper Ocean Processes Group (UOP) launched by the Woods Hole Oceanographic Institution (WHOI) focuses on the study of physical processes in the upper ocean and at the air-sea interface using moored surface buoys equipped with meteorological and oceanographic sensors. These observations enable the first accurate quantification of the annual cycle of net air-sea heat exchange and wind stress from a Southern Ocean location [34]. Twenty-three sites for the UOP were collected. The Tropical Atmosphere Ocean/Triangle Trans-Ocean Buoy Network (TAO/TRITON) [24] in the tropical Pacific, The Pilot Research Moored Array in the Tropical Atlantic (PIRATA) [35] in the tropical Atlantic, and The Research Moored Array for African-Asian-Australian Monsoon (RAMA) [23] in the tropical Indian Ocean, were the components of the Global Tropical Moored Buoy Array (GTMBA) program deployed in the tropical oceans [36]. The data provided by the three networks have undergone extensive quality control examinations to ensure that they meet stringent accuracy standards [37,38], and they have been used in studies for understanding, monitoring and forecasting ENSO and ENSO-like events, meridional modes, and monsoon variability [23]. Thirty-five sites from the three networks were selected in this study (TAO: 21, PIRATA: 7, and RAMA: 7). The OceanSITES network consists of buoys funded by individual oceanographic investigators from many different nations. The aim of the OceanSITES network is to measure, deliver and promote the use of high-quality multidisciplinary data from long-term, high-frequency observations at fixed locations in the open ocean [12]. Eight sites from OceanSITES were used. Observations from 66 moored buoy sites have been widely applied in previous studies, such as air-sea interaction studies [13,[39][40][41], ENSO prediction [42], and monsoon variability studies [23].
In general, the instruments for measuring the radiative components R ↓ g and R ↓ l were the Eppley Laboratory precision spectral pyranometer (PSP) and the precision infrared radiometer (PIR) with calibration accuracies of 2% and 1%, respectively [38]. According to previous study [12], the radiative flux measurements provided by the moored buoys are generally thought to be the highest-quality compared with those obtained from vessels; hence, the moored radiative measurements could be taken as the "ground truth" for tuning satellite measurements and assessing uncertainties in the surface flux estimates from satellites and models. For meteorological factors, RH and T a were usually measured by the Rotronic Instrument corporation MP-101 temperature-humidity probe with accuracies of 2.7% and 0.2 K, respectively [38]. These radiative and meteorological sensors were deployed at a height of 3-5 m above the sea surface at the moored buoy site. To be different, the SST was measured under approximately 1 m of the ocean surface using Sea Bird Electronics (SBE37/39). As Table 2 shows, the observations from the different observation networks have different characteristics resulting from the different measurement instruments, data-processing procedures, quality-control policies, and so on; hence, all radiative and meteorological measurements were preprocessed first in the present study, and the details are given below.

Radiative Measurements
As mentioned above, only two downward radiations (R ↓ g and R ↓ l ) were routinely measured at the ocean surface. Hence, the ocean surface R n at moored buoys could be calculated by the following formulas [43,44]: where α is the daily ocean surface shortwave broadband albedo, ε ocean is the daily ocean surface broadband emissivity, and σ is the Stefan-Boltzmann's constant (5.67 × 10 −8 W·m −2 ·K −4 ). From Equations (4) and (5), R ↓ g , R ↓ l , and SST can be measured directly, while the other two key parameters, α and ε ocean needed to be carefully defined for R n calculation. In most previous studies, α and ε ocean were usually highly simplified in parameterizations using specific constants [28,45]. However, several recent studies have pointed out that the variations in α and ε ocean cannot be ignored, especially in the ocean surface energy balance calculation [12,46,47]. In this study, the α dataset from Feng et al. [48] during 2000-2018 and the ε ocean (8-13.5 µm) dataset from Cheng, Cheng, Liang, Niclòs, Nie and Liu [46] during 1981-2018, both with a 10 km spatial resolution at a daily scale, were used. The α dataset was generated using satellite reflectance data via MODIS and other ancillary information via MERRA reanalysis data based on a three-component reflectance model of ocean water, and its validation accuracy was generally consistent and satisfactory with the previous parameterization scheme [48]. The ε ocean dataset was generated based on its relationship with wind speed, and its validation accuracy was less than 0.003 under wind-free conditions at a daily scale, which was lower than the measurement uncertainty of 0.004 [46]. Hence, the employment of the two datasets was more reasonable. Note that the multiyear average daily α from 2000 to 2018 was used as the value for α before 2000 by considering the very small annual changes in α [49]. Furthermore, the buoy SST observations only represent the bulk SST because of the radiative cooling of the ocean surface skin; therefore, it is necessary to adjust the bulk SST measurements into ocean surface skin SST before R n calculation. The bulk SST was closely related to the skin SST [50], the relationship was affected by wind speed, and a constant correction offset can be used as long as wind speed exceeds 6 m·s −1 [51]. In fact, 83% of all samples were in a condition with wind speeds exceeding 4 m·s −1 , and the moored, measured bulk SST was directly adjusted to the skin SST by using a cooling bias of 0.17 K, as Vanhellemont [52] suggested. In this study, to obtain the daily ocean surface R n at moored buoys, all measurements (R ↓ g , R ↓ l , and adjusted SST) were aggregated into daily means from hourly values without any missing data and then fed into Equations (4) and (5). For those measurements at higher observation frequencies (i.e., 10 min), they were aggregated into hourly means first only when the available observations were more than five in one hour. Note that only the measurements flagged as high-quality were used in this study. Finally, the in situ ocean surface R ↓ g , R ↓ l , and R n were obtained for modeling and validation afterwards. However, it is difficult to collect all inputs required in Equations (4) and (5), for instance, most moored buoys usually provide only one of the downward shortwave or longwave radiation measurements; therefore, an effective but more practical model for ocean surface R n estimation is needed to meet the requirements of various studies and applications.

Meteorological and Other Ancillary Information
As Table 1 shows, three meteorological variables (T a , RH, and ASTD) were also needed in this study. It is known that T a is an essential parameter for radiation estimation. However, the difference between T a and SST (ASTD) cannot be ignored for ocean surface net longwave radiation estimation. According to Fung et al. [53], the omission of ASTD could result in large net longwave radiation estimation errors in summer and winter for midlatitude and subarctic regions. Hence, the ASTD was taken into account for ocean surface R n modeling in this study, as suggested by Clark, Eber, Laurs, Renner and Saur [45] and Bunker [54]. In addition, RH was also considered in this study because it is usually used as an indicator to capture the contribution of longwave radiation from moisture in the lower atmosphere [55][56][57][58]. The T a and RH could be obtained directly from the moored buoy sites, and the ASTD was computed as the difference between the skin SST and T a . As radiative measurements, the daily T a , RH, and ASTD were also calculated from their hourly values.
In addition to the radiative and meteorological variables, the other two variables, Dd and CI (see Table 1), were needed for ocean surface R n estimation model development. Dd is defined as the difference between sunrise and sunset, and the sunrise and sunset times were computed by using the National Oceanic and Atmospheric Administration's (NOAA) Solar Calculation algorithm (NOAA Solar Calculator, www.esrl.noaa.gov/gmd/grad/solcalc/ (accessed on 15 October 2021)).
CI is the ratio of the surface downward solar radiation to the extraterrestrial solar radiation (DSR toa ), and it is usually used to represent sky conditions and atmospheric turbidity, which have significant influences on the surface radiative budget balance [20,41,59]. The CI value ranges from 0 to 1, and a larger CI indicates a clearer sky. CI has been successfully employed in previous studies [1,[60][61][62] and is utilized in this study. CI can be calculated as following [1]: d r = 1 + 0.033 cos 2πDOY 365 (8) where G sc is the solar constant (1366.7 Wm −2 ), d r is the inverse relative distance from the Earth to the Sun, ω s is the sunset hour angle (rad), ϕ is the latitude (rad), δ is the sun declination (rad), and DOY is the day of the year. The CI calculated using the daily R ↓ g measurements at moored buoys was used for modeling.

Uncertainty Analysis
According to a previous study [63], the errors in the measurements at the ocean surface should be taken into account in modeling due to their relatively significant influences. Moreover, the in situ ocean surface R n cannot be directly obtained at the moored buoy sites; therefore, its uncertainty should also be considered for better reference. According to the error propagation rules [64], the standard error of calculated R n can be derived as: where S is the standard error, Y is the dependent variable, X j is the independent variable in the formula, and N is the number of X. For example, if Y is the ocean surface R n calculated by Equations (4) and (5), then N is equal to 5. This method is a general technique for estimating the moments of the variables under consideration. The first-order, secondmoment, uncertainty analysis was used because R n is a simple algebraic sum of radiative components. An implicit assumption of Equation (11) is that errors in different variables are independent. In this study, Equation (11) was applied to estimate the uncertainty of the in situ calculated ocean surface R n , and it was extended and revised to estimate the uncertainty of the newly developed model by excluding the errors of all inputs.

CERES SYN1deg_Ed4A Product
For comparison, the daily R ↓ g and R ↓ l from the CERES SYN1deg_Ed4A product was utilized in the newly developed R n estimation models. The CERES experiment provides Earth radiation budget data, including the top-of-atmosphere (TOA) radiation flux and surface radiation flux (https://ceres.larc.nasa.gov (accessed on 15 October 2021)). The CERES SYN1deg series product was calculated using the Langley Fu-Liou radiative transfer model from the input geostationary (GEO) radiance, MODIS and GEO cloud properties, Goddard Modeling and Assimilation Office (GMAO) atmospheric profiles, and MODIS aerosols [14]. Daily R ↓ g was computed by averaging the hourly R ↓ g and extracted according to the locations of the moored buoys.

A New Ocean Surface R n Empirical Estimation Scheme
As described above, the R n at the land surface is closely related to the R ↓ g or R ↓ l determined by daytime duration [1,60,65], which has not been fully explored at the ocean surface. Hence, the relationships between the ocean surface R n and its two components (R ↓ g and R ↓ l ) were analyzed first before modeling. Based on the results, empirical models for R n estimation were built and validated afterwards. In this study, the model accuracy is indicated by three statistical indices: root-mean-square-error (RMSE), mean bias error (bias) and R-squared (R 2 ).
3.1. Relationships between the Ocean Surface R n and R ↓ g /R ↓ l Similar to the study by Chen, He, Jiang and Liang [65], the variations in the relationships between the ocean surface R n and R ↓ g /R ↓ l with the length of daytime duration are explored. For convenience, a normalized index named the length ratio of daytime (LRD, 0~1), calculated using the ratio of the total daytime hours to the entire day (24 h), was used to represent daytime duration. The LRD was divided into 20 bins (by 0.05 increments), and then, the Pearson correlation (R) between the in situ daily ocean surface R n and R ↓ g /R ↓ l was calculated for each LRD bin and is shown in Figure 2.  The results indicate that the ocean surface R n was highly correlated with the R ↓ g (R ↓ l ) when the LRD was greater (less) than 0.4. Hence, it is more reasonable to estimate the ocean surface R n under different conditions defined by the LRD.

Conditional Model Development
Many empirical models have been successfully developed for R n estimation at the land surface directly from the R ↓ g [1,60,66] when the two have a high correlation. Therefore, two conditional models were developed in this study in which the ocean surface R n was estimated mainly from the ocean surface R ↓ g when the LRD was greater than 0.4 (hereinafter referred to as the Case 1 model), and mainly from the ocean surface R ↓ l when the LRD was equal to or less than 0.4 (hereinafter referred to as the Case 2 model). Table 3 shows the number of samples used for the two conditional model developments. Note that the samples for Case 2 modeling were much less than those for Case 1, because most moored buoys were located in low-and mid-latitude areas, where the LRD is usually larger than 0.4.  Figure 3 shows the scatter plot between the daily ocean surface R n and R ↓ g when the LRD was greater than 0.4. The daily ocean surface R n was highly linearly related to R ↓ g , with an R 2 value of 0.925, RMSE of 15.501 Wm −2 , and bias of 0 Wm −2 (Equation (12)). However, the linear relationship was not as good when the ocean surface R n value was smaller than 100 Wm −2 . R n = 0.846R ↓ g − 27.305 (12) To improve the accuracy of the estimated R n , other variables (see Table 1) were introduced into Equation (12). Therefore, the relationships between the model residuals ∆R n1 (the difference between the original and model estimates) from Equation (12) and these parameters were explored one by one. After multiple experiments, it was found that the RH, CI, and ASTD were related to the variations in ∆R n1 in different formats, as shown in Figure 4.
The scatter plots between the ∆R n1 and the CI, RH, and ASTD are shown in Figure 4a,c,e, respectively. To better display their relationships, the corresponding box plots are given in Figure 4b,d,f, respectively, in which the mean of the ∆R n1 and its standard error of the mean (SEM) were calculated for each bin of the CI, RH, and ASTD (by 10% increments), respectively. Specifically, the ∆R n1 varied with the CI with a quadratic relationship (Figure 4b), and with the RH (Figure 4d) and ASTD (Figure 4f) in an approximate linear relationship. Hence, the CI, RH, and ASTD were introduced into Equation (12) in their proper formats as Equation (13): The training results of Equation (13) are shown in Figure 5. Compared with the fitting results from Equation (12) (Figure 3), the accuracy of the R n estimations from the new model was significantly improved as the R 2 value increased from 0.925 to 0.971 and the RMSE decreased from 15.501 Wm −2 to 9.650 Wm −2 . Therefore, Equation (13) was determined as the Case 1 model. Figure 6 gives the scatter plot between the in situ daily ocean surface R n and R ↓ l when the LRD ≤ 0.4. The daily ocean surface R n is also linearly related to R ↓ l , although not as closely as the relationship between R ↓ g and R n in Case 1. Similar to the Case 1 model development, a simple linear regression model (Equation (14)) taking R ↓ l as the only input, was built first. The fitting accuracy of Equation (14) was accepted with an R 2 value of 0.558 and RMSE of 12.815 Wm −2 .   As in previous work, the relationships between the model residuals ∆R n2 from Equation (14) and the parameters in Table 1 were explored one by one, and two parameters (CI and Dd) were found to be linearly related to the variations in the ∆R n2 , as shown in Figure 7. Therefore, the CI and Dd were directly introduced into Equation (14) as Equation (15): Figure 8 shows that the fitting accuracy of Equation (15) was significantly improved compared with the results from Equation (14), by increasing the R 2 value from 0.558 to 0.868 and decreasing the RMSE from 12.815 to 7.014 Wm −2 , which indicated that the CI and Dd were very helpful in ocean surface R n estimation when the LRD was small. It is known that the upward longwave radiation is mainly determined by SST; hence, SST should also be taken into account for estimating R n under the Case 2 condition. Furthermore, a close linear relationship was found between the model residuals from Equation (15) (∆R n3 ) and the SST, as Figure 9 shows. Therefore, the final Case 2 model is given in Equation (16):  The fitting results of Equation (16) are shown in Figure 10. Compared with Figures 6 and 8, the results from Equation (16) were closer to the in situ ocean surface R n , with an R 2 value of 0.984 and an RMSE of 2.428 Wm −2 . Therefore, Equation (16) was determined as the Case 2 model.

Case 2 model
Overall, two conditional models (Case 1 and Case 2) for ocean surface R n estimation were successfully built based on the measurements collected from more than 60 globally distributed moored buoys and in the formats of Equation (13) and Equation (16). Specifically, when the length of daytime duration is sufficiently large (LRD > 0.4), the daily ocean surface R n is mainly determined by R ↓ g , and RH, CI and ASTD also make contributions. For other cases (LRD ≤ 0.4), the Case 2 model was formulated as a linear function of R ↓ l , Dd, CI and SST, and the fitting accuracies of both cases were satisfactory.

Validation Results
The overall accuracy of the ocean surface R n estimations derived from the Case 1 and Case 2 models was validated against all independent validation samples (see Table 2), and the results are shown in Figure 11. Based on the results, the new ocean surface R n empirical estimation scheme performed very well, with an R 2 value of 0.972, an RMSE of 9.768 Wm −2 , and a bias of −0.092 Wm −2 . Moreover, the performances of the Case 1 and Case 2 models were also satisfactory, with R 2 values of 0.972 and 0.979, RMSEs of 9.805 and 2.824 Wm −2 , and biases of −0.095 and 0.346 Wm −2 , respectively. For better illustration, validation was conducted at each moored buoy site. Figure 12 shows the spatial distribution of the validation accuracy for each moored buoy, represented by the RMSE value. Generally, the two models had better performances over the oceans in low-latitude areas than in mid-and high-latitude areas, especially over the open ocean. It was also found that the two conditional models performed more poorly over the ocean with high SSTs, where the variations in wind and atmospheric convection are more active because of intensive ocean-atmosphere interactions [67].  To further explore the influences of latitude and distance to the coastline on the performances of the newly developed models, the variations in the validation accuracy with the two factors were analyzed, and the results are presented in Figures 13 and 14. By considering the different sizes of the samples, the relative RMSE (rRMSE) was used. In order to better illustrate the influence of the coastal seas on the model performance, the results in Figure 13 were calculated based on measurements only from the buoys over open seas, while the one in Figure 14 based on all available measurements. Figure 13 shows the performances of the new models in areas with different latitudes, which were classified into nine groups by latitude, in 10 degree intervals from 40 • S to 50 • N. The rRMSE and magnitude of the bias were both much lower over the oceans in low-latitude areas than in mid-and high-latitude areas. Similarly, the performances of the new models were examined by the distance to the coastline, from 0 to 300 km, as shown in Figure 14, and the closer to the coastline, the worse the validation accuracy of the new models. The results in Figures 13 and 14 are consistent with the conclusions derived from Figure 12. On the one hand, this was possibly caused by the much higher number of samples available in low-latitude areas and open ocean for modeling, because the performance of the empirical model may tend to reflect highly sampled conditions. On the other hand, the atmospheric conditions, such as cloud characteristics, water vapor content, and aerosol types, are different for the areas near the coast and in low-latitude regions, which has not been thoroughly considered in the newly developed models. Specifically, the atmosphere over coastal ocean areas is usually affected by continental aerosols transported by favorable winds, which introduce more enhanced aerosol direct-radiative effects than over the open ocean [68,69].   Therefore, the direct validation results demonstrated the overall satisfactory performance of the newly developed models for ocean surface R n estimation. However, the accuracy of the estimated ocean surface R n obtained from the new models over the areas near the coast or at mid-high latitudes could be decreased because of the very limited number of samples available for modeling in these regions.

Model Sensitivity Analysis
To determine the effects of the input parameters on the estimated R n from the Case 1 and Case 2 models, a global sensitivity analysis was conducted by utilizing SimLab software (http://simlab.jrc.ec.europa.eu (accessed on 15 October 2021)). All inputs into the Case 1 model (R ↓ g , RH, CI and ASTD) and Case 2 model (R ↓ l , CI, Dd and SST) were entered into SimLab separately, and then 4000 combinations of the inputs for the two models (2000 in Case 1 and 2000 in Case 2) were generated by assuming a normal distribution between the lower and upper bounds for each variable. Then, 4000 ocean surface R n estimations were obtained from the Case 1 and Case 2 models with these generated combinations of inputs. Afterwards, all combinations of inputs and the corresponding ocean surface R n estimations were used for sensitivity analysis using the Fourier amplitude sensitivity test (FAST) method [70] on the SimLab platform. The total sensitivity index (TSI), which indicated the total contribution of a factor to the output variance when the interactions of other factors were also computed, was employed to determine the sensitivity of each factor. Table 4 lists the TSI of each input in the Case 1 and Case 2 models. For the Case 1 model, R ↓ g played the most important role, with the largest TSI, 86.9%, and the second most sensitive variable was the CI with a TSI of 5.5%, while for the Case 2 model, R ↓ l was the most sensitive variable with a TSI of 64.3%, and the second most sensitive variable was also the CI with a contribution of 17.7% of the variance in the outputs. Therefore, the performances of the two newly developed models depended highly on the accuracies of the ocean surface R ↓ g and R ↓ l measurements, especially the R ↓ g , because of the importance of the CI in the two models. To better understand the importance of the accuracy of R ↓ g and R ↓ l in the two conditional models, the validation accuracy of the estimated ocean surface R n from the newly developed conditional models using the daily R ↓ g and R ↓ l from CERES SYN1deg_Ed4A (CERES R ↓ g and CERES R ↓ l ) were also examined. According to previous studies, the radiation product from CERES SYN1deg_Ed4A has been acknowledged to be one of the most reliable products for ocean surfaces [71]. In this study, the accuracy of CERES R ↓ g under the Case 1 condition was also validated against all independent validation samples with an R 2 value of 0.818, an RMSE of 27.467 Wm −2 , and a bias of 0.971 Wm −2 . Additionally, the accuracy of CERES R ↓ l under the Case 2 condition was validated with an R 2 value of 0.63, an RMSE of 15.672 Wm −2 , and a bias of 0.066 Wm −2 . Afterwards, the validation results of the estimated daily ocean surface R n with the CERES R ↓ g and CERES R ↓ l (named CERES_based R n for the Case 1 model and CERES_based R n for the Case 2 model) are shown in Figure 16a,b, respectively. Note that the validation samples (14,185) for Case 1 were less than those in Figure 11 (15,952) because CERES_R ↓ g has only been available since 2000. Hence, the validation accuracy of the Case 1 model with the in situ R ↓ g measurements since 2000 was calculated with an R 2 value of 0.974, an RMSE of 9.035 Wm −2 , and a bias of −0.955 Wm −2 . Overall, the performances of the Case 1 and Case 2 models both worsened noticeably after replacing the in situ radiation measurements with those from CERES R ↓ g and CERES R ↓ l . Specifically, for the Case 1 model, the accuracy of the daily ocean surface R n estimates worsened considerably using CERES R ↓ g , with an R 2 value decreasing from 0.974 to 0.808 and RMSE increasing from 9.035 to 24.461 Wm −2 , respectively, while for the Case 2 model, the validated R 2 value decreased from 0.979 to 0.453 and RMSE and bias increased from 2.824 to 14.825 and 0.346 to 1.113 Wm −2 , respectively. Therefore, reliable and accurate daily ocean surface R ↓ g and R ↓ l datasets are necessary for generating an accurate daily ocean surface R n dataset.

Error Analysis
At the ocean surface, the errors in the moored buoy measurements mostly resulted from buoy rocking motions, sensor tilting due to wind and currents, solar contamination, thermal gradients in the dome and case temperatures, and dome contamination at sea (e.g., salt spray crystallization, bird guano, etc.) [37,72,73], and these errors in measurements have not been taken into account in the common validation (Section 4.1). Hence, the performance of the newly developed empirical models should be assessed by eliminating the influences of measurement errors in the inputs. A quantitative analysis of the error propagation of the new models was conducted in this section. In addition, the uncertainty in the in situ ocean surface R n was estimated for comparison.
First, the variables related to error analysis are defined in Table 5.
By combining with the calculation of RMSE, we obtained ∆ P-O~( 0, 9.805) in the Case 1 model and ∆ P-O~( 0.346, 2.824) in the Case 2 model. Therefore, the S m could be estimated if the uncertainty of S p − S O (S (S P −S O ) ) was known. According to the error propagation calculation (Equation (11) in Section 2.1.4) and the information provided in Table 6, S (S P −S O ) was 2.70 Wm −2 and 6.84 Wm −2 in the Case 1 and Case 2 models, respectively. Therefore, S m in the Case 1 or Case 2 model can be calculated as Therefore, the "true" uncertainty of the R n model was 10.17 Wm −2 and 6.84 Wm −2 . For comparison, the uncertainty of the calculated in situ R n using Equations (4) and (5) was estimated as 6.23 and 2.55 Wm −2 under the Case 1 and Case 2 conditions, respectively. Hence, the true uncertainty of the new model was comparable to the calculated in situ R n , but takes more accessible buoy measurements as inputs, which demonstrated that the performance of the new model was satisfactory.

Conclusions
Ocean surface R n is one of the most essential parameters in determining the ocean heat flux, which plays a vital role in controlling the stability of the Earth's climate system. However, R n , and four other radiative components are not routinely measured at the ocean surface; hence, a simple but effective ocean surface R n estimation method is necessary for practical use. In this study, a new LRD-based empirical scheme for daily ocean surface R n estimation was successfully developed based on downward shortwave or longwave radiation and meteorological measurements collected from more than 60 globally moored buoys.
The LRD-based empirical scheme includes two conditional models named Case 1 (R n~f (R ↓ g , CI, RH, ASTD)) and Case 2 (R n~f (R ↓ l , CI, Dd, SST)). Specifically, when the length of daytime duration was longer than 9.6 hours (LRD > 0.4), the ocean surface R n was primarily determined by R ↓ g , while the ocean surface R n was primarily determined by R ↓ l for the other cases. After validation against the in situ ocean surface R n calculations, the overall accuracy of the new scheme was very good with an R 2 value of 0.972, an RMSE of 9.768 Wm −2 , and a bias of −0.092 Wm −2 , and the performance of the two conditional models was also satisfactory, with R 2 values of 0.972 and 0.979, RMSE values of 9.805 and 2.824 Wm −2 , and bias values of −0.095 and 0.346 Wm −2 , respectively. Afterwards, the "true" performance of the Case 1 and Case 2 models was estimated by eliminating the measurement errors in the inputs, which were also very good, with uncertainties of 10.17 and 6.84 Wm −2 , respectively. These were comparable to the uncertainties of the calculated in situ ocean surface R n of 6.23 and 2.55 Wm −2 under Case 1 and Case 2 conditions, respectively, but with a simpler format and more easily accessible inputs. Therefore, the new empirical LRD-based scheme for ocean surface R n has strong potential to be widely used in the near future.
However, some drawbacks in the newly developed models should be noted. For example, the two models performed worse in areas near the coastline and at high latitudes because of the limited number of samples collected within these regions, and the accuracy in the R n estimates strongly depends on the accuracies of R ↓ g and R ↓ l . Furthermore, the models only work as long as the CI is available, which indicates that the new models cannot work when the daytime length is short (LRD = 0.3 in this study). More efforts should be made to improve the performance and robustness of the newly developed models in future research.