After the upgrade of its meteorological radiation observation systems in 2004, CERN had all its sensors meet the WMO standards, whereby high-precision and stable long-term data could be obtained. In order to gain long-term, high spatial resolution PAR data, the "hybrid model" was used to calculate historical broadband solar radiation based on CMA routine observations (such as sunshine duration, pressure, temperature and relative humidity), MODIS AOD data and NASA/GSFC O_{3} data. Data observed from the 39 CERN sites were used to develop the model for calculating PAR under all-sky conditions. Then the estimation model was combined with clearness index, solar elevation angle and sunshine duration to obtain the historical daily PAR data of the 724 routine stations of CMA from 1961 to 2014. Figure 2 shows the data acquisition process.

**Figure 2
**Data acquisition processThe "hybrid model" developed by Yang Kun et al.^{16–17} is a semi-physical and semi-empirical model, which retains the simplicity of the Ångström model while containing the physical process of radiation transmission. The process can be described as bellow. Due to Rayleigh scattering, ozone absorption, aerosol absorption and scattering, water vapor absorption and permanent gas absorption, solar radiation would be attenuated when it goes through the atmosphere in clear sky. The five transmittance functions are expressed as \({\tau }_{r}, {\tau }_{oz}, {\tau }_{a}, {\tau }_{w}, {\tau }_{g}\), respectively. Solar beam radiative transmittance (\({\tau }_{b}\)) and solar diffuse radiative transmittance (\({\tau }_{d}\)) can be calculated by Equations (1) and (2):

\[{\tau }_{b}={\tau }_{r}{\tau }_{oz}{\tau }_{a}{\tau }_{w}{\tau }_{g}\]

(1)

\[{\tau }_{d}=0.5\text{ }\left[{\tau }_{oz}{\tau }_{g}{\tau }_{w}\left(1-{\tau }_{a}{\tau }_{r}\right)\right]\]

(2)

The five types of transmittance can be calculated by ground pressure, atmospheric precipitation, concentration of ozone volume and atmospheric turbidity β. Detailed calculation process can be found in Yang et al.^{16–17} Daily surface radiation (\({R}_{clear}\)) under clear sky can be obtained by Equation (3):

\[{R}_{clear}={\int }_{{t}_{1}}^{{t}_{2}}\left({\tau }_{b}+{\tau }_{d}\right){R}_{0}dt\]

(3)

where \({R}_{0}\) is the solar radiation at the top of the atmosphere, \({t}_{1}\) is sunrise and \({t}_{2}\) is sunset.

The cloud transmittance (\({\tau }_{c}\)) can be gained by the sunshine duration, and the parameterization scheme is shown in Equation (4), where *n* is actual sunshine duration, \({N }_{s}\) is the maximum possible sunshine duration, the length of time for which solar direct normal irradiance exceeds a threshold value of 120 W·m^{-2} in clear-sky conditions. Under cloudy skies, daily surface radiation \({R}_{s}\) can be calculated by Equation (5).

\[{\tau }_{c}=0.2505+1.1468n/{N}_{s}-0.3974{\left(n/{N}_{s}\right)}^{2}\]

(4)

\[{R}_{s}={\tau }_{c}\cdot {R}_{clear}\]

(5)

Because PAR variation is consistent with the change of solar radiation, some scholars have established a PAR parameter estimation model with clear sky index (\({K}_{s}\)) and sine of solar elevation angle (\(\mu\)). Clear sky index can be obtained by Equation (6).

\[{K}_{s}={R}_{s}/{R}_{0}\]

(6)

The relationship between hourly mean PAR and \(\mu\) is expressed by the power law equation (7), where \({PAR}_{m}\) is the maximum value of PAR per \({K}_{s}\), and *e* determines how PAR varies with \(\mu\). The dependence of \({PAR}_{m}\) on \({K}_{s}\) is described by Equation (8), where a, b, c, d are the fitting parameters.

\[PAR=PA{R}_{m}×{\mu }^{e}\]

(7)

\[PA{R}_{m}=a+b×{K}_{s}+c×{K}_{s}^{2}+d×{K}_{s}^{3}\]

(8)

The daily value of PAR is calculated by Equation (9). \(\overline{{K}_{s}}\)is the ratio of the daily surface radiation to the daily extraterrestrial solar irradiance; \(\overline{\mu}\) is the daily average value of the sine of the solar elevation angle from sunrise to sunset; \({t}_{d}\) is the daily sunshine duration; A, B, C, D, E are parameters relying on climatic zones.

\[PA{R}_{daily}=\left(A+B×\overline{{K}_{s}}+C×{\overline{{K}_{s}}}^{2}+D×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{E}×{t}_{d}\]

(9)

The parameterization schemes for PAR data reconstruction in China's eight regions are shown in Table 1.

**Table 1
**PAR estimation models for the eight regions**Typical station** | **Region** | **Estimation equations** |

Hailun | NEC | \(\left(0.28+9.01×\overline{{K}_{s}}+2.03×{\overline{{K}_{s}}}^{2}-1.89×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.19}×{t}_{d}\) |

Beijing Forrest Station | NCP | \(\left(0.03+10.57×\overline{{K}_{s}}-4.44×{\overline{{K}_{s}}}^{2}+3.37×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.06}×{t}_{d}\) |

Shapotou | NC | \(\left(0.24+10.18×\overline{{K}_{s}}+1.43×{\overline{{K}_{s}}}^{2}-1.78×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.24}×{t}_{d}\) |

Fukang | NWC | \(\left(0.44+7.97×\overline{{K}_{s}}+5.84×{\overline{{K}_{s}}}^{2}-5.42×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.12}×{t}_{d}\) |

Lhsasa | TP | \(\left(2.67-5.83×\overline{{K}_{s}}+30.42×{\overline{{K}_{s}}}^{2}-19.37×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.14}×{t}_{d}\) |

Yanting | SWC | \(\left(0.20+9.22×\overline{{K}_{s}}+1.34×{\overline{{K}_{s}}}^{2}-1.43×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.25}×{t}_{d}\) |

Dinghu Mountain | SEC | \(\left(0.07+9.47×\overline{{K}_{s}}-2.10×{\overline{{K}_{s}}}^{2}+2.26×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.06}×{t}_{d}\) |

Lake Dong | EC | \(\left(0.18+9.26×\overline{{K}_{s}}+0.91×{\overline{{K}_{s}}}^{2}-1.01×{\overline{{K}_{s}}}^{3}\right)×{\overline{\mu }}^{1.18}×{t}_{d}\) |