A Semi-Empirical All-Sky Photosynthetically Active Radiation Based on Sentinel-2 for High-Resolution Land Surface Analysis

1. Introduction

Solar radiation is the fundamental energy source that sustains terrestrial life and, within it, photosynthetically active radiation (PAR) defines the spectral window through which plants capture energy for photosynthesis. By directly regulating carbon assimilation, vegetation growth, and ecosystem productivity, PAR plays a central role in plant functioning from individual crop canopies to the global terrestrial biosphere [1]. Its magnitude and temporal variability influence not only the efficiency with which vegetation converts light into biomass, but also the capacity of ecosystems to respond to seasonal dynamics, climatic variability, and environmental stress [2,3,4]. Therefore, accurate characterization of PAR is essential for ecological monitoring, agricultural applications, and land surface models that seek to represent vegetation productivity in space and time [5,6].

Despite its importance, PAR remains difficult to characterize continuously across space and time. Direct measurements are limited to relatively sparse sensor networks, whereas satellite- and reanalysis-based radiation products involve trade-offs among spatial resolution, temporal frequency, latency, and physical complexity [7]. Satellite estimation of PAR has long been recognized as essential for extending observations beyond point measurements, yet most operational products are designed primarily for regional to global applications rather than fine-scale land analysis [5,8]. This limitation has become increasingly relevant as land-surface studies rely more heavily on high-resolution Earth observation data to resolve vegetation structure, ecosystem functioning, and land cover heterogeneity.

Several existing products already provide valuable PAR or radiation inputs for large-scale applications. Satellite-based products such as MODIS MCD18 [9,10] and the VIIRS VNP18/VJ118 family [11] deliver PAR at kilometer-scale resolution with sub-daily to daily temporal coverage, while global composite versions are distributed at coarser 0.05∘ grids. CERES SYN1deg products provide another widely used source of surface radiative information, including daily and monthly PAR fields at 1∘ spatial resolution [12]. NASA GMAO reanalysis products, including MERRA-2, also provide meteorological and radiative forcing fields that are widely used in productivity modeling frameworks such as MOD17 products where incident PAR is commonly approximated from incoming shortwave radiation using a fixed conversion ratio [13,14]. Reanalysis datasets such as ERA5 and ERA5-Land offer another important source of radiative forcing, providing hourly surface solar radiation fields at approximately 0.25∘ and 0.1∘ resolution, respectively, from which PAR can be approximated [15,16]. These products are highly valuable for regional to global studies, but their spatial support remains substantially broader than the scales at which field observations, agricultural parcels, and heterogeneous land-cover mosaics are often analyzed. Although incoming radiation under clear-sky conditions is generally smoother in space than land-surface properties, such coarse radiation fields may still be poorly aligned with the finer-resolution (e.g.10–30 m) vegetation variables, routinely derived from optical satellite observations.

The growing availability of Sentinel-2-derived land-surface variables [17,18,19] has increased the demand for radiative inputs that are spatially consistent with fine-resolution optical observations. While Sentinel-2 does not directly detect PAR, it provides global, systematic, high-resolution observations together with atmospheric and scene-level information that can support rapid approximation of clear-sky radiative conditions. In particular, solar geometry, aerosol optical thickness, water vapor content, and scene classification information available in Sentinel-2 Level-2A products provide the necessary inputs for a physically informed PAR derivation framework that is directly aligned with Sentinel-2 Multispectral Imager (MSI) observations. Such an approach is especially relevant for near-real-time applications, where simplicity, reproducibility, and compatibility with existing Sentinel-2 workflows are essential. In this context, the objective of a higher-resolution PAR product is not necessarily to resolve atmospheric variability at the scale of individual meters, nor to replace dedicated radiative-transfer-based radiation products. Rather, the aim is to provide PAR inputs that are spatially aligned with Sentinel-2 observations and sufficiently accurate for high-resolution land-surface analysis under clear-sky conditions. This positioning is particularly important for applications that combine PAR with fine-resolution vegetation variables, where mismatches in spatial support can introduce additional uncertainty into downstream analyses.

In this study, we present a semi-empirical framework for deriving daily clear-sky and all-sky PAR based on Sentinel-2 Level-2A imagery. The approach starts by combining analytically derived extraterrestrial radiation based on solar geometry with simplified atmospheric transmittance terms parameterized using Sentinel-2 aerosol, water-vapor, and clear-sky screening information. Then, we extend the clear-sky formulation to all-sky conditions through a cloud-transmission factor derived from cloud probability. Our estimates are validated using flux-tower observations obtained from AmeriFlux sites (spanning North and South America) and compared with well-established PAR products from satellite and reanalysis sources (i.e.MODIS MCD18C2.062 and CERES SYN1de), enabling an evaluation of the sensitivity of the framework to key modeling choices. Our main goal is to assess whether Sentinel-2 can support an accurate, rapid, high-resolution, and operationally useful PAR product for land-surface applications.

2. Materials and Methods

2.1. Semi-Empirical Clear-Sky and All-Sky PAR Estimation Framework

Figure 1 summarizes the overall workflow used to derive daily clear-sky and all-sky PAR based on Sentinel-2 Level-2A imagery and to evaluate the resulting estimates against flux-tower observations and existing reference products. The workflow begins with the acquisition of Sentinel-2 atmospheric and scene classification variables, followed by the calculation of daily extraterrestrial radiation from solar geometry. The clear-sky PAR is then estimated by applying simplified atmospheric transmittance terms and converting the resulting daily surface radiation to PAR. The clear-sky estimate is then extended to all-sky conditions using a cloud-transmission factor derived from Sentinel-2 cloud probability, which attenuates PAR according to the estimated cloudiness of each pixel. Pixel-level estimates are aggregated over tower-centered windows, 250 m radius around towers, to account for footprint representativeness and are subsequently compared with in situ PAR and satellite-based PAR products.

2.1.1. Solar Geometry and Daily Extraterrestrial Irradiance

The incoming solar radiation over a horizontal surface ($G_{oh}$) is expressed in (${\mathrm{Wm}}^{-2}$) as,

$$G_{oh}=G_{on}cos{\theta }_z$$

(1)

where $G_{on}$ is the extraterrestrial solar irradiance corrected for the time of year, and ${\theta }_z$ is the solar zenith angle. Extraterrestrial irradiance is the incoming solar radiation at the top of the atmosphere and is defined as

$$G_{on}=G_{sc}{E}_0$$

(2)

where $G_{sc}$ is the average amount of solar energy received per second by a square meter of area oriented towards normal to the Sun, also known as the solar constant, and $E_0$ is the eccentricity correction factor. While $G_{sc}$ is not a physical constant, it is approximated as 1367 ${\mathrm{Wm}}^{-2}$. The extraterrestrial irradiance $G_{on}$, on the other hand, varies seasonally primarily due to Earth’s elliptical orbit around the Sun and fluctuations in solar activity.

The eccentricity correction $E_0$, or the distance correction, is based on the inverse-square law governing radiative flux, and is represented as

$$E_0={\left({\displaystyle \frac{r_0}{r}}\right)}^2$$

(3)

where $r_0$ is the mean distance between Earth and Sun, which represents 1 Astronomical unit (1 AU), and r is the instantaneous distance. As Earth’s elliptical orbit causes the Earth-Sun distance to vary periodically over the year, this correction can be expressed as a harmonic series. Spencer [20] derived a Fourier series representation of Eq.3 by fitting harmonics to the Keplerian distance variation.

$$E_0=1.00011+0.034221cos\beta +0.00128sin\beta +0.000719cos\left(2\beta \right)+0.000077sin\left(2\beta \right),$$

(4)

where $\beta$ is the day-angle calculated by multiplying the mean orbital angular velocity with the number of day,

$$\beta ={\displaystyle \frac{2\pi (N-1)}{365}}.$$

(5)

In here, the day of the year N is 1 for January 1 and 365 for December 31 (366, if leap year).

The cosine of the solar zenith angle is formulated from spherical trigonometry as,

$$cos{\theta }_z=sin\varphi sin\delta +cos\varphi cos\delta cos\mathrm{\Omega },$$

(6)

where $\varphi$ is the geodetic latitude, $\delta$ is the solar declination, and $\mathrm{\Omega }$ is the hour angle.

The solar declination $\delta$ in degrees is approximated by

$$\delta =23.45sin\left({\displaystyle \frac{360}{365}}(284+N)\right).$$

(7)

The hour angle $\mathrm{\Omega }$ is defined as

$$\mathrm{\Omega }=15(t_{\mathrm{solar}}-12),$$

(8)

where $t_{\mathrm{solar}}$ is the solar time in hours. Solar noon corresponds to $\mathrm{\Omega }=0$.

Combining the above relationships, the instantaneous extraterrestrial irradiance over a horizontal surface becomes,

$$G_{oh}=G_{sc}{E}_0\left[sin\varphi sin\delta +cos\varphi cos\delta cos\mathrm{\Omega }\right].$$

(9)

The total daily extraterrestrial radiation $H_0$ (J  ${\mathrm{m}}^{-2}$ ${\mathrm{day}}^{-1}$) at the top-of-the-atmosphere is given by

$$H_0={\int }_{sunrise}^{sunset}{G}_{oh}\left(t\right)dt,$$

(10)

which can be approximated by

$$H_0={\displaystyle \frac{86400}{\pi }}{G}_{on}\left[{\mathrm{\Omega }}_{ss}sin\varphi sin\delta +cos\varphi cos\delta sin{\mathrm{\Omega }}_{ss}\right],$$

(11)

where ${\mathrm{\Omega }}_{ss}$ is the sunset hour angle (radians).

2.1.2. Atmospheric Attenuation

Clear-sky surface radiation can be expressed as the TOA extraterrestrial irradiance attenuated by the atmosphere, which means atmospheric radiative transfer between the TOA and the ground must be modeled. Consequently, the clear-sky surface radiation can be formulated as,

$$H=H_0{t}_{\mathrm{atm}},$$

(12)

where $t_{\mathrm{atm}}$ accounts for the cumulative effects of atmospheric transmittance. A common approximation is to represent the cumulative atmospheric transmission as the product of transmittance terms associated with major attenuation processes, including Rayleigh (molecular) scattering ($t_R$), aerosols ($t_a$), and water vapor absorption ($t_w$).

$$t_{\mathrm{atm}}=t_R{t}_a{t}_w.$$

(13)

These components can significantly impact the transmission of radiation through the atmosphere, which affects the amount that ultimately reaches the surface. We adopt a computationally efficient parameterization for Rayleigh transmittance derived from the clear-sky model of Bird and Hulstrom [21]:

$$t_R=exp\left(-0.0903·m^{0.84}\right),$$

(14)

where m is the optical air mass. In its simplest form, optical air mass can be approximated from the solar zenith angle as,

$$m\left({\theta }_z\right)={\displaystyle \frac{1}{cos{\theta }_z}}.$$

(15)

However, this geometric approximation becomes less accurate at large solar zenith angles because it does not account for the curvature of the atmospheric path. To mitigate this limitation, optical air mass was estimated using the Kasten–Young formulation [22],

$$m_r\left({\theta }_z\right)={\left[cos{\theta }_z+0.50572{\left(96.07995-{\theta }_z\right)}^{-1.6364}\right]}^{-1},$$

(16)

where ${\theta }_z$ is the solar zenith angle in degrees and $m_r$ is the relative optical air mass.

Because atmospheric attenuation depends on air pressure, relative air mass was converted to pressure-corrected air mass as

$$m=m_r{\displaystyle \frac{P}{P_0}},$$

(17)

where $P_0$ is the standard sea-level pressure and P is the local surface pressure. Surface pressure was estimated from elevation using an exponential height correction,

$$P=P_{0}exp\left(-{\displaystyle \frac{z}{H}}\right),$$

(18)

where z is elevation and H is the atmospheric scale height.

Aerosol effects are represented using a Beer-Lambert transmittance based on the aerosol optical depth at 550 nm (${\mathrm{AOD}}_{550}$) and the pressure-corrected air mass $m\left({\theta }_z\right)$:

$$t_a=exp\left(-{\mathrm{AOD}}_{550}m\right).$$

(19)

This formulation provides a simplified, computationally efficient estimate of aerosol-induced attenuation, while neglecting the spectral variability of aerosol extinction across the PAR range and the partitioning of surface irradiance into direct and diffuse components. Although PAR spans 400–700 nm, we approximate aerosol effects using ${\mathrm{AOD}}_{550}$ as a practical broadband proxy. The wavelength 550 nm is commonly treated as representative because it lies near the center of the visible window and, by extension, the PAR domain. Consequently, many widely used satellite and reanalysis aerosol products report AOD primarily at 550 nm as a standard reference variable (e.g. MODIS-based retrievals and global reanalyses/assimilations such as MERRA-2 and CAMS) [23,24].

In addition to Rayleigh scattering and aerosol attenuation, atmospheric water vapor contributes to the reduction of incoming radiation reaching the surface. Considering a simplified clear-sky parameterization based on Bird and Hulstrom [21], water-vapour transmittance was represented as

$$t_w=1-{\displaystyle \frac{0.077{\left(\mathrm{TCWV}m\right)}^{0.38}}{1+0.077{\left(\mathrm{TCWV}m\right)}^{0.38}}},$$

(20)

where TCWV is the total column water vapor. This term provides a computationally efficient approximation of water-vapour-induced attenuation within the clear-sky formulation, while remaining consistent with the simplified treatment adopted for the other atmospheric transmittance components.

2.1.3. Daily PAR Derivation

Using the daily extraterrestrial radiation term $H_0$ defined in Section 2.1.1, clear-sky surface PAR was approximated as 48% of incoming shortwave radiation [5], such that,

$${\mathrm{PAR}}_{\mathrm{cs}}=0.48H_0{t}_{\mathrm{atm}},$$

(21)

where $t_{\mathrm{atm}}$ is the effective atmospheric transmittance under clear-sky conditions.

To derive daily PAR, atmospheric attenuation was not represented solely by the noon solar geometry. Instead, optical air mass was evaluated at multiple representative solar zenith angles during the daylight period, and the resulting transmittance estimates were combined to approximate an effective daily atmospheric transmission. This approach provides a computationally efficient representation of diurnal variation in the atmospheric optical path.

To extend the formulation from clear-sky to all-sky conditions, a cloud-transmission factor ($t_{\mathrm{cloud}}$) was introduced based on Sentinel-2 cloud probability:

$${\mathrm{PAR}}_{\mathrm{all}}={\mathrm{PAR}}_{\mathrm{cs}}{t}_{\mathrm{cloud}},$$

(22)

where

$$t_{\mathrm{cloud}}=1-f_{\mathrm{cloud}}\left(1-{\tau }_{\mathrm{cloud}}\right).$$

(23)

Here, $f_{\mathrm{cloud}}$ is the cloud fraction, probability-weighted cloudiness factor derived from Sentinel-2 cloud probability layer. Cloud probabilities below 20% were set to zero to suppress low-confidence cloud detections, while probabilities equal to or above this threshold were scaled from percent to fractional units:

$$f_{\mathrm{cloud}}=\left\{\begin{array}{cc}0,\hfill & P_{\mathrm{cloud}}<20,\hfill \\ P_{\mathrm{cloud}}/100,\hfill & P_{\mathrm{cloud}}\ge 20,\hfill \end{array}\right.$$

(24)

where $P_{\mathrm{cloud}}$ is the Sentinel-2 cloud probability in percentage. The parameter ${\tau }_{\mathrm{cloud}}$ is the cloud transmittance, defined as the fraction of PAR that still reaches the surface through cloud cover. In this study, ${\tau }_{\mathrm{cloud}}$ was assigned an empirical value of 0.5. Hence, the clear-sky estimate represents potential PAR under cloud-free conditions, whereas the all-sky estimate approximates the actual PAR received at the surface after accounting for cloud attenuation.

2.2. Sentinel-2

The Copernicus Sentinel-2 mission is a sun-synchronous, near-polar orbiting, multispectral satellites designed to provide high-resolution, high-revisit optical imagery for monitoring land surfaces and coastal/inland waters. Sentinel-2 Level-2A products deliver atmospherically corrected surface reflectance (Bottom-Of-Atmosphere, BOA) derived from Level-1C top-of-atmosphere measurements, together with scene classification (SCL) and quality information. As part of the Level-2A processing chain, auxiliary atmospheric fields are produced and used to model atmospheric absorption and scattering effects. Among these auxiliary layers, the aerosol optical depth (AOD) characterizes the column aerosol loading that influences path radiance and visibility, while the total column water vapor represents the amount of atmospheric moisture that drives strong absorption features (notably in the near-infrared). Providing AOD and water vapour auxiliary datasets alongside surface reflectance improves traceability of the correction, supports uncertainty-aware interpretation, and can be valuable for downstream applications such as vegetation monitoring, aquatic remote sensing, and atmospheric screening.

2.3. Other Satellite-Based PAR Products

Other satellite-based PAR products considered in this study included MCD18C2.062 and CERES SYN1deg. MCD18C2.062 is a MODIS Terra–Aqua combined Level-3 global all-sky photosynthetically active radiation (PAR) product distributed on the 0.05^o Climate Modeling Grid. It provides incident PAR in the 400–700 nm waveband as eight 3-hourly layers per day (GMT 00:00, 03:00, 06:00, 09:00, 12:00, 15:00, 18:00, and 21:00) [10]. In this study, daily PAR was derived by multiplying each 3-hourly PAR layer by the 3-hour interval length and summing the resulting eight interval totals for each day.

We also considered the CERES SYN1deg product family, which provides radiative fluxes on a coarser 1^o grid and explicitly separates direct and diffuse components of surface radiation. CERES SYN1deg combines hourly CERES and geostationary top-of-atmosphere (TOA) fluxes, MODIS/VIIRS and GEO cloud properties, MODIS/VIIRS aerosol information, and Fu–Liou radiative transfer calculations to derive surface and in atmosphere fluxes consistent with the CERES-observed TOA fluxes. The SYN1deg product is available in multiple temporal aggregations from hourly products to monthly [12,25]. We considered three CERES-based PAR quantities: surface all-sky PAR, surface clear-sky PAR, and TOA all-sky PAR. Daily CERES PAR values for these products were calculated by summing the direct and diffuse PAR components provided by the product.

2.4. Fluxnet Towers

To evaluate the satellite-derived PAR estimates, we used eddy-covariance flux towers from the AmeriFlux network (Figure 2). In addition to providing continuous measurements of ecosystem carbon exchange, particularly net ecosystem exchange (NEE) and respiration-related fluxes, these sites also collect key meteorological variables, including incoming shortwave radiation and incoming photosynthetic photon flux density (PPFD), which are directly relevant for PAR assessment. The daily data used in this study span the period from 1 April 2017, corresponding to the start of the Sentinel-2 Level-2 data record, to 31 December 2024. From the more than 506 AmeriFlux stations screened in this study, 173 sites were retained after filtering for the availability of PPFD measurement, which was required to derive an in situ PAR reference. Stations lacking this variable or falling outside the study temporal coverage were excluded from the analysis.

Daily in situ PAR was derived from AmeriFlux PPFD_IN measurements. The PPFD_IN variable is reported in the AmeriFlux daily product as the daily mean photosynthetic photon flux density ($\mu$mol photons ${\mathrm{m}}^{-2}$ ${\mathrm{s}}^{-1}$), aggregated from half-hourly observations. Daily mean PAR in W ${\mathrm{m}}^{-2}$ was obtained using the conversion factor of 4.57 $\mu$mol photons ${\mathrm{J}}^{-1}$ [26].

3. Results

3.1. Performance of the PAR Derivation Framework Against AmeriFlux

To assess how each component of the proposed framework improves PAR estimation, we evaluated four model formulations of increasing complexity against AmeriFlux observations. The first experiment represents TOA PAR, in which PAR is derived solely from extraterrestrial radiation without atmospheric correction. The second experiment introduces atmospheric attenuation to estimate clear-sky surface PAR, while adopting a constant air-mass approximation. The third experiment further refines the clear-sky formulation by allowing air mass to vary with solar zenith angle, thereby representing changes in optical path length more realistically during the day. Finally, the fourth experiment extends the framework to all-sky surface PAR by incorporating a cloud-transmission factor derived from cloud probability. This stepwise design makes it possible to isolate the contribution of each modeling component and to evaluate how the transition from TOA PAR to clear-sky and all-sky surface PAR affects agreement with tower observations.

Throughout all stages of the experiment, daily PAR estimates were generated from Sentinel-2 overpasses coincident with the flux tower locations. As Sentinel-2 acquisitions are temporally irregular and do not provide continuous daily coverage, tower-based PAR observations were matched to the dates of valid Sentinel-2-derived PAR estimates.

3.1.1. TOA PAR from Solar Geometry

The purely geometric formulation, which estimates PAR directly from extraterrestrial radiation without accounting for atmospheric attenuation, represents TOA PAR rather than surface PAR. As expected, this first-order approximation substantially overestimated the in situ PAR observations from AmeriFlux towers across all seasons (Figure 3), because it does not include radiative losses associated with Rayleigh scattering, aerosols, water vapour, or clouds. In all observations, the TOA-based formulation showed a positive bias of $+63.80W{m}^{-2}$, an RMSE of $70.58W{m}^{-2}$, and a correlation of $r=0.866$. Consequently, the TOA-based solution should be interpreted as a potential upper-bound estimate of incident PAR under idealized conditions, rather than as a realistic representation of the radiation actually reaching the land surface.

Despite this strong positive bias, the solar-geometry-based formulation was still able to reproduce the broad seasonal pattern of incoming radiation, indicating that astronomical controls such as solar declination, day length, and solar elevation account for a substantial fraction of the temporal variability in PAR. This is reflected in the relatively strong overall correlation with in situ measurements, even though the magnitude of PAR was consistently overestimated.

Furthermore, the seasonal performances of TOA PAR (Figure 3b-e) show that the magnitude of overestimation was not uniform throughout the year. The estimation errors were the smallest in winter, with a bias of $+37.75W{m}^{-2}$ and $r=0.862$, but increased markedly during spring ($+67.80W{m}^{-2}$, $r=0.782$) and especially summer ($+87.28W{m}^{-2}$, $r=0.551$), when the TOA formulation diverged most strongly from the surface observations. In autumn, the bias decreased again to $+48.63W{m}^{-2}$ and the correlation recovered to $r=0.856$. In summer, the spread around the one-to-one line also increased, which indicates weaker agreement with in situ PAR.

3.1.2. Clear-Sky Surface PAR with Atmospheric Attenuation

Introducing atmospheric attenuation under a constant-air-mass assumption substantially improved the PAR estimates relative to the purely geometric TOA formulation, yielding a clear-sky surface PAR product that more closely matched AmeriFlux observations (Figure 4). In this experiment, optical air mass was treated as constant over the day and calculated from the solar zenith angle at solar noon. By accounting for radiative losses associated with Rayleigh scattering, aerosols, and water vapour, the systematic overestimation was strongly reduced. Overall, the clear-sky surface formulation showed a bias of $+14.28W{m}^{-2}$, an RMSE of $25.93W{m}^{-2}$, and a correlation of $r=0.877$.

Compared with the TOA case, the scattering of estimations shifted much closer to the one-to-one line, confirming that atmospheric attenuation accounts for a substantial share of the gap between extraterrestrial and surface PAR. However, the remaining positive bias indicates that the clear-sky formulation still tended to overestimate PAR, particularly under conditions where a constant air-mass approximation is too restrictive.

Seasonally, the best agreement was obtained in winter, with a bias of $+7.98W{m}^{-2}$, an RMSE of $15.99W{m}^{-2}$, and $r=0.869$. Spring and autumn remained reasonably consistent, with biases of $+16.29W{m}^{-2}$ and $+8.33W{m}^{-2}$, respectively, whereas summer, again, showed the largest difference from tower observations, with a bias of $+20.99W{m}^{-2}$ and a reduced correlation of $r=0.605$. Thus, although atmospheric attenuation substantially improved the surface PAR estimates, the constant-air-mass clear-sky formulation remained less reliable during high-radiation periods.

3.1.3. Clear-Sky Surface PAR with Varying Air Mass

Allowing optical air mass to vary with solar zenith angle further improved the clear-sky surface PAR estimates compared with the constant-air-mass formulation (Figure 5). This step reduced the overall positive bias from $+14.28W{m}^{-2}$ to $+6.84W{m}^{-2}$, while also lowering the RMSE to $23.05W{m}^{-2}$. The overall correlation remained similar ($r=0.873$), which indicates that the main benefit of the varying-air-mass formulation was a better representation of PAR magnitude rather than a major change in temporal agreement. The improvement was strongest in winter and autumn, where the remaining biases were close to zero ($+0.41W{m}^{-2}$ and $+0.69W{m}^{-2}$, respectively). These seasons also showed relatively low errors, with RMSE values of $13.89W{m}^{-2}$ in winter and $16.78W{m}^{-2}$ in autumn. In spring, the formulation still showed a moderate positive bias of $+9.19W{m}^{-2}$, while summer remained the most challenging season, with the largest bias ($+13.53W{m}^{-2}$), highest RMSE ($28.91W{m}^{-2}$), and lowest correlation ($r=0.591$).

The varying-air-mass formulation provides a more realistic clear-sky surface PAR estimate by accounting for changes in atmospheric optical path length during the day. However, the remaining positive bias, particularly in spring and summer, indicates that clear-sky atmospheric attenuation alone is not sufficient to represent actual surface PAR under real conditions. Therefore, this limitation motivates the so called all-sky formulation, where cloud attenuation is explicitly included.

3.1.4. All-Sky Surface PAR with Cloud Attenuation

Introducing cloud attenuation substantially reduced the systematic bias in the Sentinel-2 PAR estimates (Figure 6). Across all seasons, the all-sky formulation produced a near-zero mean bias of $-1.11W{m}^{-2}$, compared with the positive bias observed for the clear-sky formulation. The overall correlation remained high ($r=0.873$), while the RMSE was $23.54W{m}^{-2}$.

The seasonal results show that the effect of cloud attenuation was most important in spring and summer. In spring, the bias was reduced to $-0.84W{m}^{-2}$, while in summer it decreased to $+0.54W{m}^{-2}$. These values show that the all-sky correction effectively removed much of the seasonal overestimation associated with the clear-sky formulation. However, the error spread remained relatively large in both seasons, with RMSE values of $26.76W{m}^{-2}$ in spring and $29.79W{m}^{-2}$ in summer. This suggests that the cloud attenuation term improved the mean behaviour of the retrieval, but did not fully resolve day-to-day variability in cloud conditions.

Winter and autumn showed the strongest agreement with the tower observations. The winter estimates had a bias of $-2.24W{m}^{-2}$, an RMSE of $13.63W{m}^{-2}$, and a correlation of $r=0.875$. Autumn showed similar behaviour, with a bias of $-2.33W{m}^{-2}$, an RMSE of $17.55W{m}^{-2}$, and a correlation of $r=0.876$. The smaller errors in winter and autumn likely reflect, at least in part, the lower absolute PAR range during these seasons. However, this should not be interpreted as evidence that cloud-related uncertainty is necessarily weaker in winter, since cloud effects depend on local cloud regime and atmospheric conditions.

Overall, the all-sky formulation provided a more balanced representation of surface PAR than the clear-sky formulation. The main improvement was the reduction of systematic overestimation, particularly during spring and summer, when ignoring cloud attenuation led to larger positive biases. At the same time, the remaining scatter indicates that a single Sentinel-2 observation and its associated cloud probability cannot fully represent the complete daily evolution of cloud cover. Therefore, the all-sky product should be interpreted as a high-resolution approximation of daily surface PAR that improves the clear-sky estimate, while still retaining uncertainty under variable cloud conditions.

3.2. Comparison with Existing PAR Products

3.2.1. Benchmarking Against MODIS and CERES

Table 1 summarizes the agreement between the evaluated satellite-derived PAR products and the flux tower observations for all records combined and for each season. Overall, the Sentinel-2 all-sky PAR product performed better than the corresponding Sentinel-2 clear-sky formulation, mainly because the cloud adjustment reduced the systematic overestimation. Across all seasons, the mean bias decreased from $6.84W{m}^{-2}$ for the clear-sky product to $-1.11W{m}^{-2}$ for the all-sky product, while the overall correlation remained unchanged ($r=0.87$). This indicates that the cloud-correction step improved the accuracy of the Sentinel-2 retrieval, especially under conditions where clear-sky assumptions would otherwise lead to overly high PAR estimates.

Among the satellite PAR products, MODIS all-sky PAR showed the strongest overall agreement with the flux tower measurements. It had lower RMSE and MAE than both Sentinel-2 products and also achieved a higher correlation (r = 0.93). This is not surprising, since MODIS PAR is already designed as a physically based all-sky product and includes sub-daily information. Even so, the Sentinel-2 all-sky estimates performed reasonably well, especially considering that it provides much finer spatial detail than MODIS.

The comparison with CERES further placed the Sentinel-2 and MODIS results into context. CERES surface all-sky PAR at the daily scale showed very strong agreement with the tower observations, with lower overall error than the Sentinel-2 and MODIS products. In contrast, CERES surface clear-sky PAR and CERES TOA all-sky PAR showed much larger positive biases across all seasons. This was expected, as these quantities are not directly comparable to in situ surface all-sky PAR measurements: clear-sky PAR represents potential radiation in the absence of clouds, whereas TOA PAR does not account for atmospheric attenuation at the surface.

Seasonal statistics revealed clear differences in product performance. The benefit of the Sentinel-2 all-sky correction was most evident in spring and summer, when the clear-sky formulation showed substantially larger positive bias. In winter and autumn, the differences between the two Sentinel-2 products were smaller. The lowest correlations for Sentinel-2 were found in summer, likely reflecting stronger day-to-day cloud variability and the difficulty of representing daily integrated PAR from satellite observations that do not fully resolve the diurnal cycle. The MODIS all-sky PAR remained more stable throughout the seasons, with consistently lower bias and error metrics than the Sentinel-2 products.

Overall, these results show that the Sentinel-2 all-sky product is a clear improvement over the Sentinel-2 clear-sky formulation.Although the daily MODIS and CERES surface all-sky products achieved lower overall errors, the Sentinel-2 all-sky estimates provide substantially finer spatial detail and allow PAR to be derived directly from Sentinel-2 observations.

3.2.2. Comparison of Sentinel-Based PAR Estimate with Global Products

Figure 7 compares the spatial patterns of the Sentinel-2 all-sky PAR estimates with MODIS and CERES for representative scenes from each season. The comparison with CERES all sky PAR should be interpreted mainly as a coarse-scale magnitude evaluation rather than a pixel-level spatial comparison. CERES shows broad, coarser resolution patterns, and therefore cannot capture the fine scale cloud variability observed by Sentinel-2. In this case, MODIS provides a more useful intermediate comparison, as it retains some sub-regional variability while still being much coarser than Sentinel-2. The Sentinel-2 product generally follows the broad MODIS patterns in cloud-affected regions, but with sharper spatial transitions and more local detail.

In general, the Sentinel-2 product shows substantially finer spatial detail than both reference products, as expected from its higher spatial resolution. The comparison also shows that the cloud adjustment introduces spatially coherent reductions in PAR, particularly in areas where cloud features are visible in the Sentinel-2 RGB images. This pattern is most evident in spring and autumn, where lower PAR values in the Sentinel-2 product broadly coincide with cloud-affected regions and are also reflected, although more coarsely, in the MODIS PAR fields.

Seasonal differences are also visible. Winter shows relatively limited spatial contrast compared with the other seasons, while spring and autumn show stronger cloud-related spatial variability. Summer presents a different case, with large areas close to clear-sky conditions and comparatively high PAR values. In this situation, the Sentinel-2 product appears more strongly influenced by the clear-sky component of the algorithm, and the comparison suggests that remaining differences with MODIS and CERES may be larger where cloud attenuation is limited.

4. Discussion

In this study, we evaluated two Sentinel-2 PAR formulations: a clear-sky version and an all-sky version that includes a cloud adjustment. The contrast between the Sentinel-2 clear-sky and all-sky products is one of the most informative parts of the analysis. The clear-sky formulation can be understood as a potential PAR estimate under cloud-free conditions, based on solar geometry and atmospheric transmittance terms. In that sense, it is not wrong. However, it is not directly comparable to actual tower PAR under real daily conditions, because cloud attenuation is a major control on the amount of radiation that reaches the surface. This explains why the clear-sky version consistently shows positive bias, especially in brighter seasons.

The all-sky formulation improves this by introducing an empirical cloud-transmission term based on cloud probability. The clearest gain is the reduction in systematic overestimation. When cloud effects are ignored, the Sentinel-2 clear-sky product tends to produce PAR values that are too high relative to the flux tower observations. After the cloud-adjustment step is applied, this bias is reduced substantially while the overall temporal agreement remains similar. This step is necessary if Sentinel-2 is to be used as a practical PAR product rather than only as a clear-sky benchmark. The strongest benefit appears in spring and summer, when the difference between clear-sky and all-sky estimates becomes much larger. This seasonal behavior is physically sensible, since cloud effects during high-radiation periods can lead to large differences in incoming radiation. In contrast, the smaller differences in winter and autumn suggest that under lower-radiation conditions the impact of the cloud correction is still present, but less dominant.

Daily MODIS all-sky PAR and daily CERES surface all-sky PAR both perform better overall against the tower observations (see Table 1). In fact, this is not unexpected considering both products were designed specifically for surface radiation applications and benefit from more direct treatment of all-sky radiative conditions at coarse spatial resolution. Even so, the Sentinel-2 all-sky estimates remain meaningful as they provide finer spatial detail and allow PAR to be derived directly from the Sentinel-2 observation framework itself. This makes them useful in settings where PAR needs to be analyzed together with Sentinel-2 bands, vegetation indices, and other Sentinel-2-derived variables (e.g. GPP, NPP, or biomass) in a consistent way.

Although the Sentinel-2 estimates are generally within the seasonal range in comparison with MODIS and CERES, clear differences remain still. In areas that appear close to clear-sky conditions, Sentinel-2 often produces relatively high PAR values (see Figure 7). This suggests that some residual overestimation may remain when the cloud-transmittance correction is relatively less effective. This behavior is consistent with the structure of the algorithm, since low cloud probability leaves the estimated PAR largely controlled by the clear-sky radiation component. In contrast, where stronger cloud attenuation is present, the Sentinel-2 spatial patterns become more comparable to the MODIS patterns, although Sentinel-2 preserves finer local variability.

4.1. Limitations of the Method

An important limitation of the Sentinel-2 approach is that the sensor observes the surface only at one moment in the morning. This means that the cloud adjustment used here is based on atmospheric conditions near overpass time rather than on the full evolution of cloud cover throughout the day. That limitation is especially relevant because the target variable in this study is daily, not instantaneous PAR. A morning cloud scene does not necessarily represent afternoon conditions, and in many environments the cloud field can change substantially over the course of the day.

The cloud-transmission, although based on clear-sky index which is already a known factor [27], depends on an empirical term. It performs well in the current comparison, but it is not derived from a direct physical relationship between cloud probability and cloud optical transmittance. This means that its transferability outside the calibration domain remains uncertain. While AmeriFlux provides a strong testbed because of its broad latitudinal coverage, the transferability of the parameterization across all global biomes and climate regimes may still require further refinement to achieve higher accuracy at the global scale.

Furthermore, the atmospheric transmittance is still approximated rather than fully integrated over the entire daylight period. Although increasing the number of sampled hour angles improves the representation of daytime variation, the method still relies on a simplified daily treatment of the atmosphere. In particular, AOD and TCWV are derived from the conditions observed in instantaneous Sentinel-2 overpass and are assumed to represent the complete daytime atmospheric state. This introduces an additional first-order approximation as the sub-daily variability in aerosols, water vapor, and cloud-related atmospheric attenuation is not explicitly resolved. This is different from a full sub-daily radiative framework and likely contributes to residual errors.

Part of the residual overestimation under near clear-sky conditions may also be related to the fixed shortwave-to-PAR conversion factor used in the formulation. The fraction of broadband shortwave radiation represented by PAR is not constant, but varies with solar geometry and atmospheric composition. Therefore, adopting a fixed value of 0.48 may lead to overestimation in daily PAR values where the actual PAR fraction is lower [28]. Since low cloud probability leaves the estimate largely controlled by the clear-sky component, this source of uncertainty is expected to be most evident under clear or weakly attenuated conditions.

4.2. Future Work

Although the all-sky formulation improved the Sentinel-2 PAR estimation compared with the clear-sky PAR, the cloud transmittance term still requires a more realistic representation of daily atmospheric conditions rather than relying on an instantaneous snapshot to cloud coverage. Therefore, we plan to derive a daily clear sky index from CERES all-sky and clear-sky products and relate this coarse-scale attenuation factor to Sentinel-2 cloud probability and other atmospheric variables. This would allow CERES to provide the large-scale radiative constraints while Sentinel-2 brings the spatial detail needed for high-resolution PAR derivation.

Another aspect that requires further investigation is the incorporation of terrain effects into Sentinel-2-based all-sky PAR estimation. The fine-scale spatial patterns observed in the Sentinel-2 PAR maps indicate that local landscape structure can influence the derived PAR field, particularly in heterogeneous and topographically complex areas [29]. Terrain parameters such as slope, aspect, elevation, illumination geometry, cast shadows, and sky-view conditions can affect the amount of PAR received at the surface, especially where opposing slopes experience contrasting solar exposure. Explicitly accounting for these effects could improve the physical consistency of Sentinel-2 PAR estimates and help separate atmospheric attenuation from terrain-driven radiation variability [30].

Finally, we plan to examine the sensitivity of daily PAR estimates to the fixed shortwave-to-PAR conversion ratio. Although, a constant ratio provides a practical first order approximation, the fraction of broadband shortwave radiation represented by PAR can vary spatio-temporally with atmospheric clearness, water vapour, solar geometry, season, and site conditions. [31,32]. Therefore, a fixed value may introduce systematic deviations, particularly under weakly attenuated atmospheric conditions where the estimated PAR is mainly controlled by clear-sky radiation component.

5. Conclusions

In this study, we developed a Sentinel-2-based PAR framework with both clear-sky and all-sky formulations to enable PAR retrieval directly from Sentinel-2 observations, thereby supporting high-resolution applications such as productivity and biomass estimation. The semi-analytical approach combines Sentinel-2 AOD and TCWV information with solar geometry equations, while the all-sky formulation further incorporates cloud attenuation using Sentinel-2 cloud probability.

Based on the validation against daily PAR observations at AmeriFlux towers, we found that the all-sky formulation reduced the overestimation found in the clear-sky version and provided a more realistic representation of surface PAR under actual atmospheric conditions. Although MODIS and CERES surface all-sky products achieved lower overall errors, the Sentinel-2 approach offers finer spatial detail and direct compatibility with other Sentinel-2-derived variables.

Overall, the present results suggest that a Sentinel-2-based all-sky PAR product is both feasible and useful. It does not replace established coarse-resolution radiation products, but it offers a meaningful alternative when high spatial detail and direct integration with Sentinel-2 observations are important. In future work, we will investigate refining the all-sky formulation by improving the daily representation of atmospheric condition and cloud attenuation, incorporating terrain-driven irradiance effects, and evaluating the sensitivity of PAR estimates to the spatio-temporal changes in shortwave-to-PAR conversion factor.