Feasibility of Estimating Cloudy-Sky Surface Longwave Net Radiation Using Satellite-Derived Surface Shortwave Net Radiation

Abstract

Surface longwave net radiation (LWNR) is a vital component in the surface radiation budget. Major progress has been made in the estimations of clear-sky LWNR. However, the estimation of cloudy-sky LWNR remains a significant challenge. In this paper, a linear model (LM) and a multivariate adaptive regression spline (MARS) model were developed to estimate the cloudy-sky LWNR from a satellite-derived surface shortwave net radiation product. Spatially and temporally matched satellite data and ground-measured LWNR, which was collected at 24 sites from four networks, were used to build and validate the linear and MARS models. The effects of land cover, climate type, and surface elevation on the estimate of LWNR were also analyzed. The MARS model, incorporating the normalized difference vegetation index (NDVI) and surface elevation (H) as the inputs, had the best performance. The determination coefficient, BIAS, and root mean square error (RMSE) were 0.51, 0.01 W/m², and 26.10 W/m², respectively. The developed model, when combined with freely distributed Global LAnd Surface Satellite (GLASS) products, showed promise for producing surface LWNR and all-sky surface net radiation.

1. Introduction

Surface net radiation plays a fundamental role in determining atmospheric and oceanic thermal conditions and circulations, which shape the main characteristics of the earth’s climate. Surface net radiation is also a necessary input for land surface models that characterize the hydrological, ecological, and biogeochemical processes on earth [1,2,3]. Surface net radiation is the sum of shortwave net radiation and longwave net radiation (LWNR) [4], which can be expressed by the following equations:

$$R_n=R_{nS}+R_{nL},$$

(1)

$$R_{nS}=R_{Sd}-R_{Su},$$

(2)

$$R_{nL}=R_{Ld}-R_{Lu},$$

(3)

where $R_n$ is the surface net radiation, $R_{nS}$ is the shortwave net radiation, $R_{nL}$ is the LWNR, $R_{Sd}$ is the shortwave downward radiation, $R_{Su}$ is the shortwave upward radiation, $R_{Ld}$ is the longwave downward radiation, and $R_{Lu}$ is the longwave upward radiation. LWNR is a key component in the surface radiation budget. Accurate estimates of surface longwave radiation are required for calculating surface net radiation—especially at night. This, in turn, controls all of the surface energy budget components, including latent heat flux [5,6].

LWNR is the difference between the surface longwave downward radiation (SLDR) and the surface longwave upward radiation (SLUR). However, because of the limited number of observation sites, ground-based measurements are seriously restricted by the spatially sparse coverage that is available for mapping large-scale radiation. Satellite remote sensing is one of the most important methods used to estimate LWNR on regional and global scales. Additionally, it can provide land and atmospheric data with high spatial resolution and global coverage at lower costs in comparison to ground-based measurements [2,7].

During the past several decades, significant efforts have been devoted to estimating SLDR and SLUR using satellite remote sensing [7,8,9,10,11,12,13,14,15]. Significant progress has been made in estimating clear-sky SLDR and SLUR. For example, Cheng et al. [8,9] proposed new hybrid methods to estimate SLDR and SLUR from the Moderate Resolution Imaging Spectroradiometer (MODIS) data on a global scale, and they produced two years (2003 and 2013) of LWNR products. Under cloudy conditions, cloud cover is the main factor that regulates longwave radiation [16,17,18]. The optical–thermal satellite sensor cannot penetrate the clouds to detect the land surface signal. In addition, it is not easy to obtain cloud information, such as cloud-top temperature, cloud-top height, and cloud optical depth, which is unstable because of the great variations in cloud distributions. Therefore, the estimation of cloudy-sky LWNR remains a challenge.

Recently, estimating LWNR, including SLDR and SLUR, under cloudy-sky conditions has attracted increasing attention from the whole remote sensing community, and some achievements have been made. Many studies have shown that SLDR under a cloudy sky can be estimated using the cloud fraction and clear-sky SLDR [5,19,20,21,22,23,24]. However, these studies only focus on SLDR and it is impossible to obtain LWNR without pre-determined SLUR. Bisht et al. [25] estimated LWNR under cloudy skies using cloud-top temperature, cloud fraction, cloud emissivity, cloud optical thickness, and land surface temperature (LST) from the MOD06 level 2 cloud data product. However, the LST of MOD06 is obtained from various data sources, which include the MOD11_L2 product, the National Centers for Environmental Prediction (NCEP) gridded analysis, and Data Assimilation Office (DAO, a NASA data assimilation office) data with a spatial resolution of 5 km. The accuracy of the LST is not guaranteed. Therefore, the way to retrieve LWNR using the derived variables from the optical–thermal observations, like Bisht et al. [25], may have large uncertainties.

Studies have demonstrated that surface shortwave net radiation can be estimated directly from satellite data in all weather conditions [26,27,28]. Surface shortwave downward radiation is the main source of the earth’s energy. The land surface becomes warmer as it absorbs solar radiation, resulting in increasing emissions of LWNR. Therefore, surface longwave radiation has a certain correlation with shortwave radiation and can be estimated directly from shortwave radiation. Zhou et al. [29] proposed a linear model to estimate LWNR from shortwave net radiation for cloudy-sky conditions, which used the ground measurement collected from the Surface Radiation (SURFRAD) network. Furthermore, because of the influence of land cover and terrain, the normalized difference vegetation index (NDVI) and elevation were incorporated to adjust the model. The superiority of this method is that the input parameters are easily obtained with high accuracy [30]. However, it is developed with ground-based data, which limits its applications. The questions of whether this method can be adapted to satellite data on a large scale and of the level of accuracy that may be achieved remain unanswered.

The purpose of this paper is to investigate the feasibility of estimating LWNR under cloudy-sky conditions using a remote sensing surface shortwave radiation product. The rest of this paper is arranged as follows: The satellite data and ground measurements that were used are introduced in Section 2. Section 3 describes the methods, which is followed by an evaluation of the results in Section 4. Section 5 is the discussion, and the conclusions are presented in Section 6.

2. Data

2.1. Satellite Data

Global LAnd Surface Satellite (GLASS) products are a series of freely available, long-term, high-quality datasets that include leaf area index (LAI), shortwave broadband albedo (ABD), broadband emissivity (BBE), downwelling shortwave radiation (DSR), and photosynthetically active radiation (PAR) at the first stage [31]. There are twelve products in GLASS, which were gradually distributed. GLASS DSR and GLASS ABD products from 2008–2010 were used in this study to estimate surface shortwave net radiation. The GLASS DSR product, with a three-hour temporal resolution and a five-kilometer spatial resolution, was produced from multi-source remote sensing data, using a look-up table algorithm [32]. The GLASS ABD was produced by integrating two intermediate albedo products that were derived by a direct-estimation algorithm with a statistics-based temporal filter (STF) [33,34]. Thus, it was more complete and accurate than the intermediate albedo products. The spatial and temporal resolutions were 1 km and eight days, respectively [34]. The albedo was assumed to be unchanged during the 8-day period. Previous studies have noted that GLASS DSR and ABD products have higher accuracy than mainstream products [35,36]. Furthermore, the MODIS vegetation indices product MOD13A2 was also used to describe the land characteristics. MOD13A2 data were provided every 16 days at a 1 km spatial resolution as a gridded level 3 product in the sinusoidal projection. Each day during the 16-day period was set to be the same as the GLASS ABD.

2.2. Ground Measurements

The most popular instrument that is used to directly measure SLUR and SLDR is called a pyrgeometer, which has evolved into the precision infrared radiometer (PIR) [1]. It is relatively expensive and sensitive when compared to the currently used pyranometers, which are employed to measure surface shortwave radiation [20]. Therefore, the measurements of SLUR and SLDR are not as common as the measurement of shortwave radiation. To support scientific research, several surface radiation observation networks that provide continuous surface longwave radiation measurements have been established, and the ground-based data are publicly available. Ground-based measurements of SLUR and SLDR, which were collected from 24 sites in four networks, were used in this study in order to develop and validate the models. These included nine sites from the Arid and Semi-Arid Region Collaborative Observation Project (ASRCOP) network [36], six sites from the Surface Radiation (SURFRAD) network [37], six sites from the AmeriFlux network [38], and three sites from the Baseline Surface Radiation Network (BSRN) [39]. All of the data underwent a manual visual inspection before they were used.

Table 1. Description of site conditions.

Site Name	Latitude (°)	Longitude (°)	Elevation (m)	Land Cover	Time Period
Bondville 1	40.05	−88.37	213	Cropland	2008–2010
Boulder 1	40.13	−105.24	1689	Grassland	2008–2010
Fort Peck 1	48.31	−105.10	634	Grassland	2008–2010
Desert Rock 1	36.63	−116.02	1007	Desert	2008–2010
Penn State 1	40.72	−77.93	376	Cropland	2008–2010
Sioux Falls 1	43.73	−96.62	473	Cropland	2008–2010
Brookings 2	44.35	−96.84	510	Grassland	2008–2010
Canaan Valley 2	39.06	−79.42	994	Grassland	2008–2010
Fort Peck 2	48.31	−105.10	634	Grassland	2008
Morgan Monroe 2	39.32	−86.41	275	Forest	2008–2010
Wind River 2	45.82	−121.95	371	Forest	2008–2010
MissouriOzark 2	38.74	−92.20	220	Forest	2008–2010
PAY 3	46.82	6.94	491	Cultivated	2008–2010
TAT 3	36.05	140.13	25	Grassland	2008–2010
TOR 3	58.25	26.46	70	Grassland	2008–2010
Arou 4	38.04	100.46	3033	Grassland	2008–2009
Dongsu 4	44.09	113.57	970	Grassland	2008–2009
Jinzhou 4	41.15	121.20	22	Cropland	2008–2009
Miyun 4	40.63	117.32	350	Cropland	2008–2009
Naiman 4	42.93	120.70	361	Desert	2008
Tongyu grass 4	44.57	122.88	184	Grassland	2008–2009
Tongyu crop 4	44.57	122.88	184	Cropland	2008–2009
Yingke 4	38.85	100.40	1519	Cropland	2008–2009
Yuzhong 4	35.95	104.13	1965	Desert	2008–2009

¹ Surface Radiation (SURFRAD) sites; ² AmeriFlux sites; ³ Baseline Surface Radiation Network (BSRN) sites; ⁴ Arid and Semi-Arid Region Collaborative Observation Project (ASRCOP) sites.

summarizes the site information, which includes latitude, longitude, elevation, land cover, climate type, and observation period. The types of land cover at the sites include desert, cropland, grassland, and forest. The elevations of the sites vary from 25 m to 3033 m. The data period covers 2008–2010. The spatial distribution of these sites is shown in Figure 1.

3. Methods

According to Equation (1), LWNR can be expressed as follows:

$$R_{nL}=R_n-R_{ns}.$$

(4)

Previous studies have shown that surface net radiation can be estimated directly from surface shortwave radiation. Kaminsky and Dubayah explored the relationship between surface net radiation and shortwave fluxes in central Canada. They suggested that a single linear equation, using surface net shortwave radiation, can be used to estimate surface net radiation [40]. Alados et al. did similar work in a semi-arid area [41]. Therefore, the right-hand side of Equation (4) can be expressed only by $R_{ns}$. This means that we can estimate LWNR by $R_{ns}$. As surface shortwave net radiation can be precisely estimated under all weather conditions, we explored the relationship between LWNR and surface shortwave net radiation to estimate LWNR under a cloudy sky. This allows for measurements of cloud information to be bypassed, including cloud-top temperature, cloud-top height, and cloud optical depth, which is unstable. A linear model and a multivariate adaptive regression splines (MARS) model were used to predict LWNR. In addition, NDVI was also incorporated into the developed models to qualify the effect of vegetation on the estimation of LWNR.

3.1. Linear Model

The linear model (LM) is expressed as follows:

$$R_{nL}=a_1{R}_{nS}+b_1,$$

(5)

$$R_{nS}=(1-\alpha )R_{Sd}.$$

(6)

When the NDVI was considered, the linear model becomes the following:

$$R_{nL}=a_2{R}_{nS}+b_2\mathrm{NDVI}+c_2,$$

(7)

where $a_1$, $b_1$, $a_2$, $b_2$, and $c_2$ are regression coefficients; $\alpha$ is the shortwave broadband albedo derived from GLASS ABD; $R_{Sd}$ is the surface shortwave downward radiation derived from GLASS DSR; and NDVI is derived from MOD13A2.

To fit the above equations, spatial and temporal matching between satellite data and site measurements were implemented through three steps: (1) Spatial matching. The location of each site was determined by matching the site’s coordinates and the geolocation fields of the images; (2) Cloudy-sky identification. The cloudy-sky condition of the sites was determined by cloud fraction $c$. The calculation for $c$ was referenced from Carmon et al. [19]. When $c$ was greater than 0.05, the sky was considered to be cloudy; (3) Temporal matching. The timing of the satellite data and measurements was matched using the nearest-neighbor method to obtain time-matched data. Following the above processes, more than 90,000 samples were obtained. The samples were randomly divided into two parts. Two-thirds of the samples were used to fit the coefficients in Equations (5) and (7) by linear regression fitting, and one-third of the samples were used to validate the developed models.

3.2. MARS Model

Multivariate adaptive regression spline (MARS) is a form of regression analysis that was introduced by Jerome H. Friedman in 1991 [42]. MARS builds models in the following form:

$$f(x)={\displaystyle \sum _{i=1}^k{c}_i{B}_i}(x),$$

(8)

where $B_i(x)$ is a basis function, which is either a constant or has the form $\max (0,x-const)$ or $\max (0,const-x)$, and $c_i$ is a constant coefficient.

MARS is a non-parametric regression technique, in which nonlinear responses between variables are described by a series of linear segments of different slopes, each of which is fitted using a basis function. The breaks between segments are defined by a knot in a model that initially over-fits the data and then is simplified using a backwards/forwards stepwise cross-validation procedure to identify the terms that are to be retained in the final model. This produces continuous models with continuous derivatives, and has more power and flexibility to model relationships that are nearly additive or involve interactions between a few variables at most.

MARS was applied for LWNR estimation in this study. To explore the nonlinearity among LWNR, net shortwave radiation, and NDVI, we used MARS to establish the nonlinear relationship using the same samples as in Section 3.1. This was implemented on the MATLAB platform with a tool called ARESLab. First, we used the function aresev to select the number of basis functions using cross-validation. Figure 2 shows the change in mean square error (MSE) for the models of each fold as the number of basis functions increases. The ten dotted pink lines indicate the MSE for the models of each fold. The solid pink line is the mean MSE for each model size. The circle and vertical dashed lines are at the minimum of the solid lines, which represent the optimum number of basis functions. As shown in Figure 2, the number of basis functions was set to 11. We then used the two part samples from Section 3.1 to develop and validate the MARS model.

4. Results

4.1. Validation of the LM Model

Using the method described in Section 3, we obtained 92,220 samples. The coefficients in Equation (5) were fitted by linear regression using two-thirds of the samples:

$$R_{nL}=-0.12R_{nS}-11.74.$$

(9)

The BIAS and root mean square error (RMSE) were used as the primary indicators of accuracy. In addition to the BIAS and RMSE, the determination coefficient (R²) was used as an indicator to test the performance of the developed model. The training results are presented in Figure 3a. The R² was 0.34, and the BIAS and RMSE were 0.006 W/m² and 30.22 W/m², respectively. N is the number of samples for training purposes. An evaluation of the model’s accuracy was then implemented using the remaining one-third of samples that were extracted. The validation results are presented in Figure 3b. N is the number of samples for validation purposes. The R² was 0.34, and the BIAS and RMSE were 0.02 W/m² and 30.19 W/m², respectively. The LWNR values were overestimated at the low end of the LWNR range—where the ground-measured LWNR values were approximately −150 W/m²—while the LWNR values were underestimated at the high end of the LWNR range—where the ground-measured LWNR values were approximately 0. In addition, there was a ceiling for the linear model. This may have been because the surface shortwave net radiation was not enough to estimate LWNR.

In addition to clouds, many factors can affect the surface radiation budget, including near-surface air temperature and humidity, soil moisture, and land cover type. However, most of these parameters cannot be obtained directly by optical–thermal satellite under cloudy conditions. The NDVI from MOD13 was used to provide additional land surface information. We then obtained a new linear model (LM-NDVI) with two variables as follows:

$$R_{nL}=-0.12R_{nS}+28.11\mathrm{NDVI}-23.76.$$

(10)

The training results are presented in Figure 4a. The R² was 0.37, and the BIAS and RMSE were 0.013 W/m² and 29.56 W/m², respectively. One-third of the samples were then extracted and used to validate the LM-NDVI model. The evaluation results of LM-NDVI are presented in Figure 4b. The R² was 0.37, and the BIAS and RMSE were 0.08 W/m² and 29.54 W/m², respectively. Compared with the LM model, the LM-NDVI model did not improve significantly, but the distribution of scatters appears to be more reasonable. There is no obvious upper limit for LM-NDVI.

4.2. Validation of the MARS Model

The MARS model was trained with the same samples that were used in Section 4.1. Figure 5a presents the training results of the MARS model. The validation process was performed to evaluate the MARS model. These results are presented in Figure 5b. The performance of the MARS model was slightly better than that of the LM model. For the validation results, the determination coefficient was 0.36, and the BIAS and RMSE values were 0.02 W/m² and 29.86 W/m², respectively. Similar to the LM model, the obvious maximum and minimum values of the retrieved LWNR were −20 and −120, respectively.

For the same reason as in Section 4.1, NDVI was considered and inputted into the MARS model for new LWNR predictions (MARS-NDVI model). Figure 6 plots the scatter density of the ground measurements against the MARS-NDVI model-predicted LWNR using the training samples and validation samples, respectively. The validation results were almost the same as the training results. This indicates that the training samples were sufficient. The validation results indicated that the LWNR predictions were well-related to the ground measurements, with a determination coefficient of 0.46, and small BIAS and RMSE values of 0.08 W/m² and 27.49 W/m², respectively. Moreover, the MARS-NDVI model performed much better than the MARS model, indicating that the NDVI was a factor that could not be neglected and needed to be incorporated into the LWNR estimations. Meanwhile, the MARS-NDVI model had a better performance than the LM-NDVI model, which indicated that the nonlinear model was more robust than the linear model.

Vegetation has obvious seasonal variations and may have certain effects on different LWNR estimation models. Therefore, in order to further investigate the adaptabilities of different LWNR estimation models among different seasons, the samples were divided into four parts—namely, spring, summer, autumn, and winter—to test the four developed models. The results are presented in Figure 7 and summarized in

Table 2. Summary of statistics of the four developed models for spring, summer, autumn, and winter. Linear model (LM); normalized difference vegetation index (NDVI); multivariate adaptive regression spline (MARS); RMSE; BIAS.

Seasons	No. of Samples	LM		LM-NDVI		MARS		MARS-NDVI
		RMSE	BIAS	RMSE	BIAS	RMSE	BIAS	RMSE	BIAS
Spring	9275	28.05	0.16	28.04	0.17	27.25	0.20	26.30	0.28
Sumer	9287	32.37	−0.29	28.24	−0.14	31.84	−0.25	26.48	−0.22
Autumn	7194	30.31	0.45	28.94	0.33	29.62	0.49	26.86	0.43
Winter	4111	28.35	0.07	28.36	0.07	27.44	−0.05	26.45	−0.12

. The results indicated that the LWNR predictions had obvious seasonal differences, without considering the NDVI. Both the LM and MARS performed the worst in summer, with the largest RMSE of greater than 30 W/m². These models had better performances in winter. This is because less surface shortwave net radiation was transformed into longwave radiation as a result of the transpiration of the vegetation. The correlation between surface shortwave net radiation and LWNR weakened with the increase in vegetation coverage. Therefore, the seasonal differences were significantly decreased when considering NDVI.

The validation accuracy of the developed MARS-NDVI model was comparable to the accuracy of the general model from Zhou et al. [29], which is reported to have a BIAS and RMSE of −2.31 and 29.25 W/m², respectively, at nine AmeriFlux sites. The developed MARS-NDVI model is a promising method for estimating LWNR using current satellite products (e.g., GLASS SDR and ABD, MODIS, and NDVI).

5. Discussion

Earth–atmosphere interactions could be affected or modulated by many factors, such as land cover, climate type, and surface elevation. Therefore, we attempted to investigate the influence of these factors on the performances of the above four developed models in this section. First, the land surface was divided into four primary land cover types—namely, desert, cropland, grassland, and forest. Then, the accuracy of each model for the different land cover types was evaluated. The evaluation results are presented in Figure 8a,b. The statistical results are provided in

Table 3. Summary of statistics of the four developed models for the different land cover types, climate types, and surface elevations (H). Semi-arid climate (BS); desert climate (BW); temperate climates (Cf); continental climates (Df); and dry winter continental climates (Dw).

Class	No. of Samples	LM		LM-NDVI		MARS		MARS-NDVI
		RMSE	BIAS	RMSE	BIAS	RMSE	BIAS	RMSE	BIAS
Desert	1752	39.84	−0.01	39.84	−0.02	38.17	−0.13	36.92	0.17
Cropland	13,693	22.32	−0.11	22.15	−0.14	22.15	−0.09	21.55	−0.07
Grassland	14,421	29.73	0.05	29.36	0.07	29.37	0.11	27.99	0.12
Forest	873	20.62	0.01	20.56	0.00	19.60	0.02	19.27	−0.24
BS	6101	35.72	0.23	35.72	0.22	35.02	0.30	34.31	0.24
BW	1601	39.99	0.15	38.87	0.06	38.40	0.05	35.92	−0.32
Cf	4543	24.69	0.17	24.24	0.19	24.47	0.15	23.09	0.20
Df	17,710	21.06	−0.12	20.82	−0.13	20.93	−0.09	20.38	−0.07
Dw	459	45.15	0.69	33.17	0.83	43.60	0.62	30.07	−0.32
H < 500m	22,649	23.85	−0.04	23.85	−0.04	23.70	−0.01	23.03	0.01
500 < H < 1000	2990	32.77	0.28	32.71	0.26	31.93	0.36	30.33	0.33
H > 1000	5100	37.85	0.19	36.66	0.30	37.14	0.21	34.41	0.17

. All four of the models performed much worse over desert than they did over vegetated areas, with an RMSE of approximately 40 W/m². This may be related to the characteristics of high albedo and low specific heat capacity in arid areas [43].

Based on the Koppen climate classification, we divided the samples into five groups: BS (semi-arid climate), BW (desert climate), Cf (temperate climates), Df (continental climates), and Dw (dry winter continental climates). We then evaluated the performance of the four models over the different climate types. The results are presented in Figure 8c,d. All of the models performed better over humid climate types (Cf, Df) than they did over dry climate types (BS, BW, and Dw). As shown in Table 3, the RMSE was approximately 20 W/m² over humid climate types and greater than 30 W/m² over dry climate types. This was consistent with the results obtained for desert land cover.

Similar to the land cover and climate type analysis, the samples were divided into three groups according to the elevation of the sites (H < 500 m, 500 < H < 1000, and H > 1000). The accuracy of the four models at the different heights is presented in Figure 8e,f. The accuracy decreased as elevation increased (Table 3). Thus, surface elevations were incorporated into the developed MARS-NDVI model (named the MARS-NDVI-H model) to improve the model’s performance. Both the training and validation results were improved. Figure 9 shows the scatter density of the ground measurements against the MARS-NDVI-H model-predicted LWNR using the same training and validation samples that were used for the MARS-NDVI model, respectively. The determination coefficient increased to 0.51, and the BIAS and RMSE were reduced to 0.01 W/m² and 26.10 W/m², respectively. Surface elevations were also incorporated into the developed LM-NDVI model (named the LM-NDVI-H model). The performance of the LM-NDVI-H model was improved, but was poorer than the performance of the MARS-NDVI-H model.

According to the above analysis, it was more desirable to develop site-specific LWNR estimation models. However, such models may have lacked representativeness due to limited ground observations. On the other hand, the uncertainties in land cover type, climate classifications, and the adopted thresholds may have resulted in larger errors in the estimated LWNR. Therefore, we preferred to use the unified model for its simplicity and acceptable accuracy at the current stage.

Though the BIAS and RMSE of the developed MARS-NDVI-H model were acceptable, the determination coefficient was relatively low. Primarily, there were two potential error sources: (1) Cloud-base temperature and vertical structure are not considered in GLASS DSR. GLASS DSR is retrieved from optical data using a look-up table method that only considers cloud type and absorption coefficient [32]. Cloud-base temperature and vertical structure can affect the accuracy of the SLDR estimate [44]. (2) The temporal mismatch between satellite data and in situ measurements. The three-hour GLASS DSR was interpolated from the derived instantaneous DSRs using the nearest-neighbor method. The satellite data and ground measurements were also matched using the nearest-neighbor method in this study. The first error source was hard to overcome because the operational method of DSR estimation, which could incorporate cloud base temperature and vertical structure, was unavailable. As a result of an unaffordable memory requirement in the production of GLASS DSR, the retrieved instantaneous DSR was not saved. We can output the instantaneous DSR as an intermediate product in the near future, which would greatly alleviate the effect of the temporal mismatch between satellite data and in situ measurements. Therefore, the second error could be greatly reduced and the determination coefficient of the MARS-NDVI-H model could be improved.

6. Conclusions

Satellite estimates of cloudy-sky LWNR remain challenging. In this study, the feasibility of estimating LWNR using satellite-derived surface shortwave net radiation was investigated. A large number of spatially and temporally matched samples under a cloudy sky were extracted from ground-based measurements, which were collected at 24 sites from four networks and remote sensing products (GLASS and MODIS data). The extracted samples were divided randomly into two parts, in which two-thirds of the matched samples were used for model development and one-third of the matched samples were used for validation purposes. Linear and MARS models for LWNR estimation were developed and validated using in situ observations. Moreover, the adaptabilities of the above LWNR estimation models among different seasons were investigated. The influence of land cover, climate type, and elevation on these four models was also discussed.

According to the validation results that were obtained from the limited data, the LM model had the poorest performance. The determination coefficient, BIAS, and RMSE were 0.34, 0.02 W/m², and 30.19 W/m², respectively. The nonlinear MARS-NDVI-H model performed the best, with a determination coefficient, BIAS, and RMSE of 0.51, 0.01 W/m², and 26.10 W/m², respectively. This accuracy was comparable to that of a previous study that estimated LWNR using ground-measured surface shortwave net radiation. We are deploying the flux measurement sites in China to acquire a time series of SLUR and SLDR measurements, and are collecting SLUR and SLDR measurements from existing flux measurement networks. We will update the developed models when more ground measurements and instantaneous GLASS DSR become available.

Our group has produced and freely distributed high-quality GLASS DSR, GLASS ABD, and clear-sky LWNR products, and it is promising to produce the corresponding LWNR product with the developed model. Additionally, we can calculate the all-sky surface net radiation from each component of the surface radiation budget.