Assessment of BRDF Impact on VIIRS DNB from Observed Top-of-Atmosphere Reﬂectance over Dome C in Nighttime

: The Day–Night Band (DNB) imaging sensor of the Visible Infrared Imaging Radiometer Suite (VIIRS) adds nighttime monitoring capability to the Suomi National Polar-Orbiting Partnership and National Oceanic and Atmospheric Administration 20 weather satellite launched in 2011 and 2017, respectively. Nighttime visible imagery has already found diverse applications, but image quality is often unsatisfactory. In this study, variations in observed top-of-atmosphere (TOA) reﬂectance were examined in terms of nighttime bidirectional effects. The Antarctica Dome C ground site was selected due to high uniformity. First, variation of reﬂectance was characterized in terms of viewing zenith angle, lunar zenith angle, and relative lunar azimuth angle, using DNB data from 2012 to 2020 and Miller–Turner 2009 simulations. Variations in reﬂectance were observed to be strongly anisotropic, suggesting the presence of bidirectional effects. Then, based on this ﬁnding, three popular bidirectional reﬂectance distribution function (BRDF) models were evaluated for effectiveness in correcting for these effects on the nighttime images. The observed radiance of VIIRS DNB was compared with the simulated radiance respectively based on the three BRDF models under the same geometry. Compared with the RossThick-LiSparseReciprocal (RossLi) BRDF model and Hudson model, the Warren model has a higher correlation coefﬁcient (0.9899–0.9945) and a lower root-mean-square-error (0.0383–0.0487). Moreover, the RossLi BRDF model and Hudson model may have similar effects in the description of the nighttime TOA over Dome C. These ﬁndings are potentially useful to evaluate the radiometric calibration stability and consistency of nighttime satellite sensors.


Introduction
The Visible Infrared Imaging Radiometer Suite (VIIRS) is a key instrument [1] onboard the Suomi National Polar-Orbiting Partnership (S-NPP) and National Oceanic and Atmospheric Administration (NOAA) 20 (formerly Joint Polar Satellite System-1, JPSS-1), launched on 28 October 2011, and 18 November 2017, respectively. To support highquality nighttime-visible imagery, the Day-Night Band (DNB), covering wavelengths from 500 to 900 nm on VIIRS, has a specified dynamic range spanning seven orders of magnitude, making it capable of detecting Earth scenes under illumination as dim as quarter moon. This range is achieved by a three-gain-stage charge-coupled device (CCD) with four arrays of imaging detectors: the low gain stage (LGS) for daytime, the medium gain stage (MGS) for dawn and dusk, and two redundant high gain stages (HGSs) for nighttime, high gain 2 of 28 stage A (HGA) and high gain stage B (HGB). Detectors in each of the DNB arrays need to be individually calibrated [2]. The VIIRS DNB sensor is also the first of its kind to provide in-orbit radiometric calibration. However, calibration of DNB is complicated, due to its large dynamic range, ultra-high sensitivity, scan-angle-dependent field of view (FOV), and sensitivity to stray light contamination [3,4]. The calibration uncertainty of the LGS and MGS are 5-10% and 10-30%, respectively, while the uncertainty of the HGS is over 30%. To improve the quality of DNB nighttime visible imagery, uncertainty of the HGS must be reduced.
The Dome Concordia (Dome C) area (75.1 • S, 123.4 • E) on the East Antarctic plateau has been considered as an ideal ground site for calibration and validation of satellite sensors, especially for VIIRS DNB HGS. The site has high reflectance from permanent ice and snow cover, and its surface is spatially very homogeneous with a small slope and low surface roughness. The high reflectance site can reflect most moonlight from the top-of-atmosphere (TOA) and transfer it into the HGS sensors, for low-energy observation. However, while the surface stability and atmospheric conditions at Dome C site offer many advantages for DNB sensor calibration and validation, a major disadvantage is that there is a large diurnal variation in reflectance, due to bidirectional effects.
To compensate for these effects, the bidirectional reflectance distribution function (BRDF) can be used to characterize the geometrical reflecting properties of a reflecting surface. In general, BRDF models correct reflectance variations due primarily to varying incident zenith angle and viewing zenith angle (VZA). Ma et al. [5] conducted the vicarious calibration of S-NPP VIIRS DNB on the basis of the deep convective clouds (DCCs), under lunar illumination, to evaluate the radiometric calibration of DNB HGS at night. However, the bidirectional reflectance of DCC was not taken into account under the criterion of lunar zenith angle (LZA) less than 60 degrees and VZA less than 60 degrees. Furthermore, Cao et al. [6] used the Angular Distribution model to estimate the BRDF effect of the DCCs and applied the reflected lunar radiance at night, from the DCC, for the intercalibration of VIIRS DNB on S-NPP and NOAA-20. However, the applicability of the Angular Distribution model for DCC at night had not been analyzed. In the retrieval of NASA's Black Marble nighttime product suit (VNP46), Roman et al. [7] adopted the RossThick-LiSparseReciprocal (RossLi) BRDF model, to correct the surface BRDF effects caused by varying illumination conditions-namely moonlight and reflected airglow from the Earth's upper atmosphere. They thought RossLi model could capture the wide range of conditions affecting the VIIRS DNB on a global basis. Moreover, all available daytime, atmospherically corrected VIIRS DNB bidirectional reflectances over a multi-data period were collected to establish the analytical solution for the RossLi BRDF model parameter values.
However, existing BRDF models are developed primarily for daytime applications, such as analysis of vegetation [8][9][10][11], soil [9], snow [10,[12][13][14], and so on. At Antarctic Dome C site, the illumination differs between day and night. Sunlight is the main light source during the day, while sources of light at night are complex, such as moonlight, airglow, lamplight, starlight, aurora, and so on. As the natural light source at night, the moon is about 250,000 times darker than the sun [15]. Additionally, air temperatures at night at Dome C site are nearly 55 • C cooler than the daytime [16]; hence, the daytime and nighttime ground characteristics may be different due to changes in environmental factors. The popular Warren and Hudson models were derived based on in situ daytime measurements over Antarctica Dome C, and the RossLi BRDF model is a widely used kernel-driven model. The applicability of these three BRDF models to nighttime over Dome C is unknown. Moreover, the specific uncertainty is unclear under the implicit assumption that daytime models are equally applicable to nighttime ones [6,17]. Note that, in this article, two BRDF models derived based on in situ near-surface measurements of the snow surface over Antarctica Dome C are named after the model-makers' names (Warren model and Hudson model represents the model of Warren et al. [18] and Hudson et al. [19], respectively). Therefore, this study first performed a systematic analysis of angles' influence on nighttime TOA reflectance and then used this analysis to explore the applicability of surface BRDF models to nighttime TOA over Dome C, using DNB data from 2012 to 2020. This paper is organized as follows. Section 2 introduces the data, three surface BRDF models, and a TOA lunar irradiance model. Section 3 introduces the study area and the methodology to simulate nighttime radiance at satellite, using three BRDF models to explore their applicability. Section 4 presents the results and related discussion. Section 5 summarizes the work.

VIIRS DNB
VIIRS is a key instrument on-board the S-NPP and NOAA-20 which are in sun synchronous orbit of about 824 km. VIIRS is a scanning radiometer with a scanning wide of about 3044 km (the cross-track direction) and a 12 h revisit time. Its wavelength ranges from 410 to 12.500 nm. According to its resolution, VIIRS includes 5 high-resolution imagery channels (I-bands), 16 moderate resolution channels (M-bands), and one DNB. The rotating telescope assembly (RTA) mirror on VIIRS reflects radiation onto a set of CCD detectors, and one rotation of the mirror is one scan. M-bands and the DNB have 16 detectors to detect this radiation, and I-bands have 32 detectors with twice the resolution of the M-bands and the DNB [20]. According to its physical properties, VIIRS has 14 reflective solar bands (RSBs), seven thermal emissive bands (TEBs), and one DNB.
The panchromatic DNB on VIIRS has a 750 m spatial resolution, and a specified dynamic range of about 7 orders of magnitude, from 3 × 10 −9 W cm −2 sr −1 to 0.02 W cm −2 sr −1 [1]. DNB achieves this dynamic range by adopting three gain stages: LGS, MGS, and HGS. The LGS is used to observe daytime scenes. The MGS is used to observe dawn or dusk scenes near the Earth's terminator. The HGS, the average of two identical stages-HGA and HGB-is used to observe nighttime scenes. Moreover, the time delay integration (TDI) of these three gain stages has 1, 3, and 250 pixels to increase signal, respectively [21]. Basic characteristics of VIIRS and DNB are listed in Table 1. More specific details can be found in Liao et al. [1]. Furthermore, the purpose of on-orbit radiometric calibration is to determine the offset and the gain. The offset of each gain stage is determined over a dark Earth scene. The gain of LGS is derived by using the solar diffuser (SD) view, while MGS and HGS saturate when the SD is illuminated. The MGS and HGS gains are calculated based on LGS results, together with the estimated gain ratios of MGS/LGS and HGS/MGS [22]. Table 1. Basic Characteristics of Visible Infrared Imaging Radiometer Suite (VIIRS) and Day-Night Band (DNB) (data with symbol "*" and "**" cite from Liao et al. [1] and Qiu et al. [21], respectively).

Ground Measured Surface Reflectance
The multi-angle ground-measured surface reflectance of Antarctic Dome C was provided by Hudson et al. [19]. Their experiments were carried out at the top of a 32 m tower, during the austral summers of 2003-2004 and 2004-2005, which is the period of polar daytime. The measurements were made by using a Field-Spec Pro JR spectroradiometer manufactured by Analytical Spectral Devices (ASD). The ASD recorded the radiance from 350 to 2500 nm, with a 25 nm resolution. Moreover, the fiber optic input cable to the ASD was mounted in a baffle, limiting its field of view to a 15 • cone. The bidirectional reflectance measurements were made under the solar zenith angle (SZA) from 51.57 • to 87 • and the VZA from 7.5 • to 82.5 • , with a 15 • interval. The key characteristics of experiments in Hudson et al. [19] are listed in Table 2. Based on these measurements data, the Hudson model was developed by using the empirical orthogonal functions (EOFs) of the data. Moreover, these measurements were also used to calculate the coefficient of RossLi BRDF model and Warren model in this study. Table 2. Key characteristics of experiments in Hudson et al. [19].

Characteristics Specification
Location RossLi BRDF model is a semi-empirical linear kernel-driven BRDF model [23][24][25][26]. It defines reflectance as a weighted sum of an isotropic parameter and two kernels, RossThick volume scattering kernel and LiSparse Reciprocal geometric-optical scattering kernel. The volume scattering kernel is derived from radiative transfer models [27], and the geometric-optical scattering kernel is based on surface scattering and geometric shadow-casting theory [28]. Moreover, the parameters in this model have corresponding physical meaning in the simulation of scattering properties of canopies and bare soils [8]. Besides this, the RossLi BRDF model has also been found to be functional when applied to reflectance observed by MODIS and S-NPP VIIRS over the surface of Antarctic Dome C [17,29,30]. Furthermore, based on framework and relevant theories of the kernel-driven BRDF model, a snow kernel was developed to better model the anisotropic reflectance of pure snow [14]. What is more, it is also the theoretical basis of the Algorithm for Model Bidirectional Reflectance Anisotropies of the Land Surface (AMBRALS), which has been applied to produce the MODIS global albedo and BRDF product [31,32]. This model is expressed as follows: where θ, ϕ, φ, and λ represent the SZA, VZA, relative solar azimuth, and wavelength, respectively. As shown in Figure 1, the SZA and VZA are made from the z-axis. The solar azimuth angle (φ 0 ) and the viewing azimuth angle (φ v ) are measured clockwise from the north. Moreover, the relative solar azimuth (φ) is defined as the angle measured clockwise from φ 0 to φ v . K vol is the RossThick volume scattering kernel, and K geo is the LiSparse Reciprocal geometric-optical scattering kernel. The f iso , f vol , and f geo represent coefficients or weightings of isotropic scattering, volume scattering, and geometric-optical scattering, respectively.
where ωl(λ) is a weight given to each observation, and d is the degrees of freedom that is number of observations minus number of parameters, fk.

Warren Model
Reflectance was measured as 600, 660, and 900 nm from a 22 m tower, at South Pole Station, in Warren et al. [18]. The SZA varied from 67° to 89.3° over the full range of VZA, azimuth angle between Sun and view, and azimuth angle between Sun and sastrugi. Measurements were made with 15° field of view at 15° intervals in VZA and azimuth angle between Sun and view. They analyzed the variation of BRDF measurements with solar elevation angle, wavelength, sastrugi azimuth, models of sastrugi reflectance, forward peak and near-nadir views. Furthermore, an anisotropic reflectance pattern was developed by using a three-term Fourier series in azimuth angle between Sun and view: where μ0 = cosθ, μr = cosφ. Thus, a total of 12 coefficients needed are determined by a leastsquares-fitting method. With enough reflectance observations ρ(λ) made at angles (θ l , ϕ l , φ l ), the analytical solution of this model coefficient, f k , is minimization ∂e 2 /∂f k = 0 of a least-squares error function [33]: where ω l (λ) is a weight given to each observation, and d is the degrees of freedom that is number of observations minus number of parameters, f k .

Warren Model
Reflectance was measured as 600, 660, and 900 nm from a 22 m tower, at South Pole Station, in Warren et al. [18]. The SZA varied from 67 • to 89.3 • over the full range of VZA, azimuth angle between Sun and view, and azimuth angle between Sun and sastrugi. Measurements were made with 15 • field of view at 15 • intervals in VZA and azimuth angle between Sun and view. They analyzed the variation of BRDF measurements with solar elevation angle, wavelength, sastrugi azimuth, models of sastrugi reflectance, forward peak and near-nadir views. Furthermore, an anisotropic reflectance pattern was developed by using a three-term Fourier series in azimuth angle between Sun and view: where µ 0 = cosθ, µ r = cosϕ. Thus, a total of 12 coefficients needed are determined by a least-squares-fitting method.

Hudson Model
Hudson et al. [19] extended the measurements of Warren et al. [18] with a wider range of SZAs (51 • -87 • ) and a broader spectrum range, from 350 to 2400 nm, at Dome C. In Hudson et al. [8] [19]), based on different spectral region and SZAs. Each group has different equations and scaling factors ( Table 2 in Hudson et al. [19]), to perform with EOFs in Equation (9): where R is a 288 × 1 matrix, 1 is a 288 × 1 matrix of ones, U is a 288 × 2 matrix gridded at 288 angular locations, Σ is a 2×2 diagonal matrix with positive scale factors in decreasing order, and V is a 1 × 2 matrix, [v1 v2]. In this paper, the first group in Table 1 and its corresponding equations and scaling factors in Table 2 of Hudson et al. [19] were used by considering the spectral region of DNB (500-900 nm) and region of LZA.

Lunar Irradiance Model
Miller-Turner 2009 (MT2009) is a TOA lunar spectral irradiance model. It was developed for calibrating nighttime low-light measurements from VIIRS DNB, to achieve quantitative nighttime multispectral applications [34]. This model was based on the solar source observations, lunar spectral albedo data, time-varying Sun/Earth/Moon distances, and lunar phase angle (LPA). If users input a specific date and time, the model is able to output the corresponding 1 nm resolution lunar irradiance, in the range of 0.2 to 2.8 µm. Figure 2 shows five spectral irradiance curves from MT2009, changing with LPA, in 2019. All the five lunar spectral irradiances are around 14:30 in May 2019. They have a similar trend with the wavelength on the whole, but the specific irradiance value is obviously different at different LPAs.

Lunar Irradiance Model
Miller-Turner 2009 (MT2009) is a TOA lunar spectral irradiance model. It was developed for calibrating nighttime low-light measurements from VIIRS DNB, to achieve quantitative nighttime multispectral applications [34]. This model was based on the solar source observations, lunar spectral albedo data, time-varying Sun/Earth/Moon distances, and lunar phase angle (LPA). If users input a specific date and time, the model is able to output the corresponding 1 nm resolution lunar irradiance, in the range of 0.2 to 2.8 μm. Figure 2 shows five spectral irradiance curves from MT2009, changing with LPA, in 2019. All the five lunar spectral irradiances are around 14:30 in May 2019. They have a similar trend with the wavelength on the whole, but the specific irradiance value is obviously different at different LPAs.

Methodology and Study Area
In this study, first of all, the S-NPP VIIRS DNB observations and MT2009 lunar irradiance model were used to calculate the TOA nighttime reflectance of Dome C during 2012 to 2020. Second, the systematical analyses of angles' influence on the nighttime TOA reflectance were conducted, including the LZA, VZA, and relative lunar azimuth angle (RAA). Third, three BRDF models (RossLi BRDF model, Warren model, and Hudson model) are respectively applied, to explore their applicability to describe the nighttime TOA of Dome C.

Study Area
The study area is in the Antarctic Dome C site (75.1 • S, 123.4 • E, elevation 3.2 km). Figure 3 shows the location of Dome C in the southern hemisphere. This site is extremely cold (temperature down to −84 • C), is covered with uniform and flat permanent snow, is cloudless most of the time, and has very high atmospheric stability due to low water vapor content, low aerosol content, and low wind speed [16]. Thus, it has been widely used in satellite sensor calibration and validation [13,[16][17][18]29,[35][36][37][38][39][40][41][42][43][44][45][46]. What is more, Dome C is located in a high latitude area near the South Pole, so polar orbiting satellites overpass this site more frequently, and the special phenomenon of polar daytime and polar night provides a clear distinction between day and night. Moreover, there is a large diurnal variation in reflectance due to bidirectional effects [45] that is also helpful to study the impact of related angles on reflectance.

Selection of Observations
The criteria to select the S-NPP VIIRS DNB data are summarized as follows: (1) The study area is a circular area centered at 75.1 • S, 123.4 • E, with a 10 km radius.
Moreover, cloud-contaminated data are removed according to the spatial uniformity (ratio of 1 standard deviation to mean reflectance over the study area) with a threshold value of 5% [44][45][46]. (2) LPA is less than 90 • to ensure sufficient moonlight.
(3) LZA is less than 75 • to ensure sufficient moonlight, too. (4) SZA is greater than 118.4 • , to remove the influences of stray-light effects present at Dome C, during the observations [21].

Methodology and Study Area
In this study, first of all, the S-NPP VIIRS DNB observations and MT2009 lunar irradiance model were used to calculate the TOA nighttime reflectance of Dome C during 2012 to 2020. Second, the systematical analyses of angles' influence on the nighttime TOA reflectance were conducted, including the LZA, VZA, and relative lunar azimuth angle (RAA). Third, three BRDF models (RossLi BRDF model, Warren model, and Hudson model) are respectively applied, to explore their applicability to describe the nighttime TOA of Dome C.

Study Area
The study area is in the Antarctic Dome C site (75.1°S, 123.4°E , elevation 3.2 km). Figure 3 shows the location of Dome C in the southern hemisphere. This site is extremely cold (temperature down to −84 °C), is covered with uniform and flat permanent snow, is cloudless most of the time, and has very high atmospheric stability due to low water vapor content, low aerosol content, and low wind speed [16]. Thus, it has been widely used in satellite sensor calibration and validation [13,[16][17][18]29,[35][36][37][38][39][40][41][42][43][44][45][46]. What is more, Dome C is located in a high latitude area near the South Pole, so polar orbiting satellites overpass this site more frequently, and the special phenomenon of polar daytime and polar night provides a clear distinction between day and night. Moreover, there is a large diurnal variation in reflectance due to bidirectional effects [45] that is also helpful to study the impact of related angles on reflectance.

Selection of Observations
The criteria to select the S-NPP VIIRS DNB data are summarized as follows: (1) The study area is a circular area centered at 75.1°S , 123.4°E , with a 10 km radius. Moreover, cloud-contaminated data are removed according to the spatial uniformity (ratio of 1 standard deviation to mean reflectance over the study area) with a threshold value of 5% [44][45][46]. (2) LPA is less than 90° to ensure sufficient moonlight.
(3) LZA is less than 75° to ensure sufficient moonlight, too. (4) SZA is greater than 118.4°, to remove the influences of stray-light effects present at Dome C, during the observations [21]. increase of LPAs. The red "+" marks the center of the site used in this study. Furthermore, more details are listed in Table 3.

Further Data Filtering According to LPA
There are the scatter plots of our calculated nighttime TOA reflectance and LPA with uniform axes during 2012-2020 in Figure 5. It is found that the relationship is a ' There is a sharply decreasing tendency when the LPA is less than 15 • , and a relative increasing tendency when the LPA is greater than 15 • . However, the phenomenon is the least obvious in 2013, and this should be related to the SRF change of S-NPP VIIRS DNB between April and May 2013, due to mirror degradation [47]. Moreover, the points whose reflectance values are extremely high mainly gather on both sides of each graph. On the left side of each graph, where the LPA is less than 5 • or so, the extremely high reflectance should have a close relationship with the opposition effects of the moon [6,34,48]. On the right side, where the LPA is greater than 70 • or so, the signal-to-noise ratio (SNR) of VIIRS DNB HGS may decrease under low moonlight. Therefore, to cut down these bad effects, the LPA criterion was then changed to 5 • -70 • . This change significantly cut down the unexpected reflectance on both sides, especially in 2019, and then the reflectance turned to be relatively uniformly distributed with LPA. In addition, a deficiency with the MT2009 model also exits. If the lunar model can be able to account for the lunar opposition surge and lower the minimum valid value of LPA, there will be less uncertainty and more data for subsequent analyses, pointing towards an improvement.  Table 3.     Table 3. Specific characteristics of these four observations in Figure 4.

Data Processing
First, the nighttime TOA reflectance R is calculated by the following formula: where L DNB is the observation radiance from VIIRS DNB, and L MT2009 is the TOA downwelling lunar radiance and is determined by the following: where θ 0 denotes LZA, E m is TOA down-welling lunar irradiance and output from MT2009 model. The following is the calculation equation of lunar irradiance: where SRF represents the spectral response function of DNB, and I MT (λ) is the lunar irradiance spectra obtained from the MT2009 model for a specific date. Second, there are three surface BRDF models, including RossLi BRDF, Warren and Hudson model. These BRDF models were respectively combined with nighttime various angles at Dome C, including LZA, VZA, and RAA, to simulate the nighttime radiance at satellite. Then the simulated radiance results are compared with the observed values to analyze the applicability of these models at night from 2012 to 2020. Figure 6 shows the flowchart of analyzing the applicability of three BRDF models at night. Moreover, the core calculation process of the BRDF models is in the black dotted box. To be specific, the coefficients of the RossLi BRDF model and Warren model were derived by using the least squares method based on the measurements data of Dome C, and the Hudson model was originally developed based on the same measurements data. Moreover, because the measurements were made at a 25 nm resolution, the average values of model coefficients from 500 to 900 nm were used to replace the coefficients of RossLi BRDF model and Warren model in DNB. Then the simulated reflectance can be calculated under arbitrary LZA, VZA, and RAA. For the Hudson model, it can only calculate the reflectance of VZA at 15 • resolution and relative azimuth at 7.5 • resolution under a given LZA and wavelength. Therefore, at first, the wavelength was set at 1 nm resolution from 500 to 900 nm. Then the reflectance R s (λ) was interpolated over VZA and relative azimuth, using the Bicubic interpolation method. Finally, the modeled reflectance, R, of DNB is obtained by convoluting R s (λ) with S-NPP DNB SRF(λ): DNB HGS may decrease under low moonlight. Therefore, to cut down these bad effects, the LPA criterion was then changed to 5°-70°. This change significantly cut down the unexpected reflectance on both sides, especially in 2019, and then the reflectance turned to be relatively uniformly distributed with LPA. In addition, a deficiency with the MT2009 model also exits. If the lunar model can be able to account for the lunar opposition surge and lower the minimum valid value of LPA, there will be less uncertainty and more data for subsequent analyses, pointing towards an improvement. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 31 model, it can only calculate the reflectance under the VZA at 15° resolution and relative azimuth at 7.5° resolution. Because the resolutions are relatively low, the Bicubic interpolation method was used in the calculation, and this processing may add some errors to results.  In this study, the use of same measurements data theoretically unified the three BRDF models. However, there were still some different errors in the application of these models. For RossLi BRDF model and Warren model, using the average values of model coefficients from 500 to 900, nm to replace the DNB coefficients brought an error. Figure 7

The Impact of Multiple Angles on Nighttime TOA Reflectance
The angular analysis is divided into three parts: VZA, LZA, and RAA.  Figure 8 shows the distribution of nighttime TOA reflectance of the five groups in nine years (2012-2020) over Dome C in the viewing hemisphere. The radius of every circle represents VZA with 10 • increments, and the polar angle of every circle represents S-NPP VAA with 30 • increments. The first column in Figure 8(a1-a5) are the scatter diagrams which can reflect the spatial distribution of the viewing geometry and the specific value of reflectance of each observation. The second column (b1-b5) plots the average reflectance in the corresponding VZA and VAA. It is evident that VAA is always between 120 • and 330 • , so the viewing geometry is stable when S-NPP overpasses Dome C. When LZA is less than 65 • and VZA is less than 40 • , the reflectance increases with the increase of VZA. Moreover, when LZA is less than 65 • and VZA is greater than 40 • , the reflectance decreases with the increase of VZA. When LZA is greater than 65 • , the reflectance has a tendency to increase with the increase of VZA.  The five rows are at LZAs between 50° and 55° (a1,b1), between 55° and 60° (a2,b2), between 60° and 65° (a3,b3), between 65° and 70° (a4,b4), and between 70° and 75° (a5,b5), respectively. The radius of every circle represents VZA with 10° increments, and the polar angle represents S-NPP viewing azimuth angle (VAA) with 30° increments. Figure 14 shows the slope and correlation coefficient (R 2 ) of linear fitting of the previous 28 graphs Figures 9-13. These values are in the same order as the graphs in Figures 9-13, and the dotted black lines divide the results into the five groups. The following conclusions can be drawn.
(1) When LZA is less than 65°, especially when the RAA approaches 90°, the reflectance tends to decrease with the increase of VZA, which is consistent with the finding in Shao et al. [47] and Qiu et al. [49]. At 55°-60° of the LZA, the average of the slopes of linear fitting is the smallest, which means that the reflectance is relatively the most stable and little affected by the VZA. (2) When LZA is higher than 65°, the reflectance increases with the increase of VZA.
Generally, the variation is more sharply at 70°-75°, so the effect of VZA reaches the most. (3) In each group divided by LZA, the value of the slope tends to decrease at first and drops to the minimum when the RAA approaches 90°. It then increases with the in- In addition, because the reflectance varies with both zenith angles (LZA and VZA), the variable-controlling approach was used to respectively analyze them. In each group they were divided into six cases, with 30 • increments, according to the RAA (0 • -30 • , 30 • -60 • , 60 • -90 • , 90 • -120 • , 120 • -150 • , and 150 • -180 • ). In each case, the LZA and RAA were regarded as constants so as to analyze the relationship between reflectance and VZA in nine years. The results are shown in Figures 9-13, where the red line in each diagram denotes the result of linear fitting. When LZA is between 50 • and 55 • , there is only one case when RAA is in 0 • -30 • and 150 • -180 • , respectively. Thus, there are no corresponding results on that conditions. Figure 14 shows the slope and correlation coefficient (R 2 ) of linear fitting of the previous 28 graphs Figures 9-13. These values are in the same order as the graphs in Figures  9-13, and the dotted black lines divide the results into the five groups. The following conclusions can be drawn.
(1) When LZA is less than 65 • , especially when the RAA approaches 90 • , the reflectance tends to decrease with the increase of VZA, which is consistent with the finding in Shao et al. [47] and Qiu et al. [49]. At 55 • -60 • of the LZA, the average of the slopes of linear fitting is the smallest, which means that the reflectance is relatively the most stable and little affected by the VZA. (2) When LZA is higher than 65 • , the reflectance increases with the increase of VZA.
Generally, the variation is more sharply at 70 • -75 • , so the effect of VZA reaches the most.    Figure 9. The relationship between nighttime TOA reflectance and VZA, in nine years, under the selection criteria (50° < LZA < 55°, 5° < LPA < 70°, and SZA > 118.4°). The red line denotes the result of linear fitting.      Figure 13. The relationship between nighttime TOA reflectance and VZA, in nine years, under the selection criteria (70° < LZA < 75°, 5° < LPA < 70°, and SZA > 118.4°). The red line denotes the result of linear fitting. Figure 14. The slope and correlation coefficient (R 2 ) of linear fitting of the previous 28 graphs (Figures 9-13). These values are in the same order as the graphs in Figures 9-13. The dotted black lines divide the results into the five groups. Figure 14. The slope and correlation coefficient (R 2 ) of linear fitting of the previous 28 graphs (Figures 9-13). These values are in the same order as the graphs in Figures 9-13. The dotted black lines divide the results into the five groups.

Impact of LZA on Nighttime TOA Reflectance
Based on the variable-controlling approach, the nighttime TOA reflectance data from 2012 to 2020 were divided into seven groups, with 10 • increments of VZA (0 In each case, the VZA and the RAA were regarded as constants, so as to analyze the relationship between reflectance and LZA. The results are shown in Figures 15 and 16, where the red line in each diagram denotes the result of linear fitting, as well as Table 4, including the ranges of VZA and RAA, linear fitting equations, correlation coefficient (R 2 ), root-mean-square-error (RMSE), and the number of cases. When VZA is in 40 • -70 • and RAA is in 150 • -180 • , there are no corresponding reflectance results. Figure 17 shows the slope and correlation coefficient (R 2 ) of linear fitting of the previous 39 results in Figure 15, Figure 16 and Table 4, and the dotted black lines divide the results into the seven groups. The following conclusions can be drawn: (1) Only when the VZA is in 50 • -60 • and RAA is less than 30 • , and VZA is in 60 • -70 • and RAA is in 30 • -120 • , the nighttime TOA reflectance is positively correlated with the LZA. In other cases, the reflectance is negatively correlated with the LZA. (2) On the whole, as the absolute value of the slope decreases, the correlation coefficient tends to decrease.
results. Figure 17 shows the slope and correlation coefficient (R 2 ) of linear fitting of the previous 39 results in Figure 15, Figure 16 and Table 4, and the dotted black lines divide the results into the seven groups. The following conclusions can be drawn: (1) Only when the VZA is in 50°-60° and RAA is less than 30°, and VZA is in 60°-70° and RAA is in 30°-120°, the nighttime TOA reflectance is positively correlated with the LZA. In other cases, the reflectance is negatively correlated with the LZA. (2) On the whole, as the absolute value of the slope decreases, the correlation coefficient tends to decrease.     Figure 15, Figure 16 and Table 4. The dotted black lines divide the results into the seven groups.
In order to reduce the error caused by too little data, less than 10 corresponding cases were deleted. So when LZA is less than 55° and VZA is in 0°-70°, and LZA is in 70°-75° and VZA is in 60°-70°, there are no corresponding reflectance results. Figure 22 shows the correlation coefficient (R 2 ) of polynomial fitting of the previous 27 graphs Figures 18-21. These values are in the same order as the graphs in Figures 18-21, and the dotted black lines divide the results into the four groups. The following conclusions can be drawn.  In order to reduce the error caused by too little data, less than 10 corresponding cases were deleted. So when LZA is less than 55 • and VZA is in 0 • -70 • , and LZA is in 70 • -75 • and VZA is in 60 • -70 • , there are no corresponding reflectance results. Figure 22 shows the correlation coefficient (R 2 ) of polynomial fitting of the previous 27 graphs Figures 18-21. These values are in the same order as the graphs in Figures 18-21, and the dotted black lines divide the results into the four groups. The following conclusions can be drawn.

Application of Three BRDF Model at Nighttime TOA
The three angles (LZA, VZA, and RAA) have different effects on the nighttime reflectance, so it is necessary to look for a suitable nighttime TOA BRDF model of Dome C site. Figure 23 shows a comparison of the simulated nighttime radiance at satellite, respectively, using three models and the observed radiance during 2012-2020. Figure  24 shows the variation of the specific correlation coefficients and the RMSE between the simulated radiance and observed radiance over the years. The application of BRDF model at nighttime TOA can be indirectly illustrated by the correlation coefficient and RMSE. According to these results of the nine years, the following conclusions can be drawn.
(1) In Figure 23, the simulated nighttime radiance at satellite using BRDF models and the observed radiance agrees well. Almost all the results are located along the 1:1 line and show high consistency with a correlation coefficient of greater than 0.9723 and an RMSE of less than 0.0799 W·cm −2 ·sr −1 .
(2) In Figure 23, the correlation of the simulated radiance and the observed has a decreasing tendency with the increase of radiance value, especially in the results of the RossLi BRDF model and Hudson model. (3) The correlation coefficients, in descending order, each year, are Warren>Hudson>RossLi.
The RMSEs, in ascending order, each year, are Warren<Hudson<RossLi. The reason why the applicability of RossLi BRDF model is lower than the other two models may be that the accuracy of RossLi BRDF model is reduced under a large zenith angle [23]. (4) During the nine years, as shown in Figure 24, the RossLi BRDF model and Hudson model have kept a good consistency. Thus, these two models may have similar effects in the description of the nighttime TOA over Dome C.     Therefore, within the margin of error of 0.55-2.77%, these surface BRDF models are suitable for the nighttime TOA over Dome C. Warren and Hudson model both were developed exactly for Antarctic snow surface, and both consist of an ample number of parameters. In this article, based on the same ground measurements, they both have a good applicability on the nighttime TOA over Dome C with higher correlation coefficients and lower RMSEs than RossLi BRDF model. The correlation coefficients of Warren model are between 0.9899 and 0.9945 always higher than RossLi BRDF and Hudson models. Moreover, the RMSEs of the Warren model are between 0.0383 and 0.0487 W·cm −2 ·sr −1 always lower than RossLi BRDF and Hudson models.

Application of Three BRDF Model at Nighttime TOA
The three angles (LZA, VZA, and RAA) have different effects on the nighttime reflectance, so it is necessary to look for a suitable nighttime TOA BRDF model of Dome C site.

Application of Three BRDF Model at Nighttime TOA
The three angles (LZA, VZA, and RAA) have different effects on the nighttime reflectance, so it is necessary to look for a suitable nighttime TOA BRDF model of Dome C site.  (Figures 18-21). These values are in the same order as the graphs in Figures 18-21. The dotted black lines divide the results into the four groups.
in the description of the nighttime TOA over Dome C. Therefore, within the margin of error of 0.55-2.77%, these surface BRDF models are suitable for the nighttime TOA over Dome C. Warren and Hudson model both were developed exactly for Antarctic snow surface, and both consist of an ample number of parameters. In this article, based on the same ground measurements, they both have a good applicability on the nighttime TOA over Dome C with higher correlation coefficients and lower RMSEs than RossLi BRDF model. The correlation coefficients of Warren model are between 0.9899 and 0.9945 always higher than RossLi BRDF and Hudson models. Moreover, the RMSEs of the Warren model are between 0.0383 and 0.0487 W⋅cm −2 ⋅sr −1 always lower than RossLi BRDF and Hudson models.    Figure 24. Variation of the correlation coefficients and the RMSE over the years between the simulated nighttime radiance at satellite, respectively, using three surface BRDF models and the observed radiance.

Conclusions
In this study, based on the variable-controlling approach, the effect of three angles (LZA, VZA, and RAA) on nighttime TOA reflectance were analyzed over Dome C, from the year 2012 to 2020, and it was found that different angles in different ranges have different influences on the reflectance. (1) When LZA is less than 65°, the reflectance tends to decrease with the increase of VZA. When LZA is higher than 65°, the reflectance increases with VZA. (2) Only when the VZA is in 50°-60° and RAA is less than 30°, or VZA is in Figure 24. Variation of the correlation coefficients and the RMSE over the years between the simulated nighttime radiance at satellite, respectively, using three surface BRDF models and the observed radiance.

Conclusions
In this study, based on the variable-controlling approach, the effect of three angles (LZA, VZA, and RAA) on nighttime TOA reflectance were analyzed over Dome C, from the year 2012 to 2020, and it was found that different angles in different ranges have different influences on the reflectance. (1) When LZA is less than 65 • , the reflectance tends to decrease with the increase of VZA. When LZA is higher than 65 • , the reflectance increases with VZA. (2) Only when the VZA is in 50 • -60 • and RAA is less than 30 • , or VZA is in 60 • -70 • and RAA is in 90 • -120 • , the nighttime TOA reflectance is positively correlated with the LZA. In other cases, the reflectance is negatively correlated with the LZA. (3) The forward scattering and backward scattering are not symmetric at night. On the other hand, three BRDF models were, respectively, used to explore the applicability of surface BRDF models at nighttime TOA, and all of these models have an excellent and suitable performance. The correlation coefficient in descending order each year is Warren>Hudson>RossLi, and the RMSE, in ascending order, each year, is Warren< Hudson<RossLi. Especially, the correlation coefficients and RMSEs of Warren model are always higher than 0.9899 and less than 0.0487 W·cm −2 ·sr −1 , respectively.
Although the Dome C has been considered as an ideal ground site for the calibration and validation of satellite sensors, there is a large nighttime variation in reflectance, due to the impact of the BRDF. Results show that bidirectional nighttime reflectance variation is very sensitive when LZA exceeds 70 degrees or VZA is lower than 10 degrees, which is therefore not recommended to conduct the evaluation of VIIRS DNB HGS calibration across satellite over Dome C. Besides, this study provided three TOA reflectance models which can be used to track the radiometric calibration stability and consistency of nighttime satellite sensors over Dome C, especially for those scanning imagers acquiring large angular views. With comparisons between three BRDF models, the findings can be expected to improve the accuracy of radiometric inter-comparisons among various sensors in the nighttime remote sensing.