An Empirical Model of Gross Primary Productivity (GPP) and Relations between GPP and Its Driving Factors, Biogenic Volatile Organic Compounds in a Subtropical Conifer Plantation in China

Abstract

Measurements of net ecosystem exchange (NEE), solar global radiation, photosynthetically active radiation (PAR) and meteorological parameters were carried out on a subtropical conifer plantation in China from 2013 to 2016. These observations were used to develop and evaluate an empirical model of gross primary production (GPP) (EMGPP) with 3-factor and 2-factor models. Using a 3-factor model, the simulated hourly GPP values were consistent with observations with a relative bias of 9.96% and normalized mean square error values of 0.07 mg CO₂ m⁻² s⁻¹ for the scattering factor S/Q (S and Q are diffuse and global solar radiation) < 0.5 and 15.52% and 0.15 mg CO₂ m⁻² s⁻¹ for S/Q ≥ 0.5. Validations of the EMGPP for hourly, daily, monthly, and annual GPP values were carried out and showed that both 3-factor and 2-factor EMGPP models can accurately capture diurnal, seasonal and interannual variations in GPP, but most simulated GPP overestimated the observed value. When the scattering factor is not available, the 2-factor EMGPP can be used. The EMGPP using 3-factor and 2-factor models was applied to simulate GPP under all sky conditions from 2013–2016, and the estimated GPP were in reasonable agreement with the measured values and showed systematic overestimations of 31% and 29% for mean hourly GPP and 41% and 29% for annual amounts, respectively. The sensitivity test demonstrated that GPP values were more sensitive to changes in PAR than to changes in water vapor and scattering factor at low S/Q, but were more sensitive to changes in water vapor than to PAR and S/Q at high S/Q. The sensitivity test revealed some mechanisms of GPP and its related processes, including the relationships between GPP and scattering of PAR, GPP and water vapor, which were in good agreement with other observations and model studies. An empirical model based on PAR energy balance can better describe the multiple interactions between GPP and its driving factors (PAR, water vapor, S/Q). The ratio of the emissions of biogenic volatile organic compounds (BVOCs) to net ecosystem exchange clearly varied between forests in different climate zones.

1. Introduction

Global warming, especially in the Antarctic, Arctic and Tibetan Plateau during the 20th century, has attracted extensive interest and has been reported on by the Intergovernmental Panel on Climate Change (IPCC) and other studies [1,2,3,4,5]. Increases in greenhouse gas (GHG, mainly CO₂, CH₄, N₂O) emissions due to anthropogenic activities are considered a common reason [6,7]. Thus, it is important to reduce greenhouse gas emissions and accurately estimate CO₂ flux and carbon storage of terrestrial vegetation. It is known that, as a large carbon sink, the terrestrial biosphere sequesters 20–30% of global anthropogenic CO₂ emissions [8].

Two methods are used to study carbon sources and sinks. (1) Monitoring and quantifying net carbon storage by forests is necessary in the validation of carbon sequestration estimates and in assessing the possible role of these ecosystems in offsetting fossil fuel emissions. The eddy covariance (EC) technique is reliable for measuring CO₂ exchange (or NEE) over periods ranging from hours to years [9], and EC flux data are usually used in validations. (2) Models are useful tools to understand the processes and mechanisms of CO₂ flux, net ecosystem exchange (NEE) and gross primary production (GPP). Models simulate NEE and GPP and investigate interactions between NEE, GPP, and their driving factors and further potential effects on regional and global climate change [10,11,12]. Large quantities of models have been developed and used, mainly including empirical models (light use efficiency (LUE), machine learning, etc. [13,14], process models [12,15,16,17,18], and atmospheric inversion models [19,20]. Most empirical models consider limited processes, use several parameters and lack understanding of mechanisms, especially multiple interactions between vegetation and the atmosphere. Process models can consider more detailed processes in vegetation, soil and atmosphere and gain a better understanding of the processes and mechanisms than empirical models. Process-based dynamic global vegetation models (DGVMs) are also widely used to estimate natural land sinks as part of the global carbon budget (GCB). Globally, the DGVM multimodel mean estimate of global natural land sinks is consistent with the estimate of global carbon budget residual land sinks, but there is a significant spread across models [18]. In addition, these complicated models have many hypotheses and assumed parameters due to our current limited knowledge of the processes involved in carbon balance [21]. Until now, there have been large uncertainties in CO₂ flux measurements (e.g., random and system errors) and model predictions [9,22,23,24,25,26]. Therefore, internal processes and external driving factors, along with their inherent connections, should be investigated thoroughly by using another method, i.e., the principle of PAR balance over the canopy as it was used in previous empirical model studies [27]. The advantages of the PAR energy method are its ability to study multiple interactions between vegetation and the atmosphere thoroughly, its use of fewer hypotheses and variables and its shorter computing time. Despite the development of measurement techniques and models, further studies on reducing measurement and model simulation errors are still needed [28,29,30].

As a part of carbon emission into the atmosphere, biogenic volatile organic compounds (BVOCs) emitted from vegetation should be considered in global carbon balance [27,31,32,33,34,35]. They are easily oxidized to produce/contribute all kinds of gases, liquids and particles (GLPs), such as O₃, HCHO, methyl vinyl ketone (MVK) and methacrolein (MaCR), CH₃(CO)OONO₂ (PAN), secondary organic aerosol (SOA), clouds [31,36,37,38,39], using ultraviolet (UV) and visible radiation, through photochemical reactions with OH and H₂O. GLPs (gases, liquids, aerosols, clouds, etc.) attenuate UV and visible radiation and influence biological processes as well as BVOC emissions [40]. Therefore, BVOCs also play significant roles in atmosphere and climate, both as participants with high chemical reactivity and as an important bridge between the gases, liquids and particles in chemical and photochemical reactions in the atmosphere [41,42,43,44]. Thus, close attention should be paid to the bidirectional impacts of GLPs and BVOC emissions [27]. In order to calculate BVOC emissions, BVOC emission models, MEGAN (Model of Emissions of Gases and Aerosols from Nature) [36,45,46,47] and others [27,48,49,50,51] were developed from various algorithms expressing photosynthesis-based schemes and are now widely used.

China has several climate zones and diverse ecosystems. It is necessary to better understand terrestrial carbon balance and its underlying mechanisms. Many studies have been carried out in China and show large differences in the magnitude of the carbon budget [52,53,54,55]. Based on reliable flux measurements and a previously used PAR energy balance method [27], it may be feasible to develop an empirical model for simulating GPP and understanding the processes and mechanisms associated with GPP.

The aim of this study is to (1) develop an empirical model of GPP using PAR energy balance to simply the complicated CO₂ process and then grasp/express its main processes, (2) evaluate the model’s performance against the observed data from a subtropical coniferous forest in China, (3) study the responses/mechanisms of GPP to its main driving factors, (4) simulate GPP from 2013–2016 under all sky conditions to deeply understand the processes and mechanisms associated with GPP and (5) study the relationships between BVOCs and NEE in several climate zones for a full understanding of carbon balance (sources and sinks).

2. Instrumentation and Methods

2.1. Site Description

This study was carried out on a conifer plantation at Qianyanzhou Subtropical Forest Ecosystem Research Station, Chinese Academy of Sciences (Qianyanzhou Station for short) in Taihe County, Jiangxi Province in subtropical China (26°44′48″ N, 115°04′13″ E, 110.8 m). The main tree species are Pinus massoniana, Pinus elliottii, Cunninghamia lanceolata and some broadleaf trees, with coverages of 42.0%, 49.6%, 5.8% and 2.6% of this landscape, respectively. These trees were planted in 1985. The leaf area index is 3.6 m² m⁻² [56]. The average canopy height is 18 m. The mean slope of this study region is 2.8–13.5°. The shrubs are mainly Loropetalum chinense, Adinandra millettii and Lyonia compta [35]. The mean annual precipitation is 1485.1 mm, and the annual temperature is 17.9 °C [57]. The annual quantity of global solar radiation was 4579 MJ m⁻² in 2013 and 4496 MJ m⁻² in 2014, and the annual quantity of photosynthetically active radiation (PAR) was 7998 mol m⁻² in 2013 and 7656 mol m⁻² in 2014 [35]. Most of the atmosphere at the Qianyanzhou site had high GLP loads, i.e., most S/Q (S and Q are diffuse and global solar radiation, respectively) values were larger than 0.5, and their hourly averages from May 2013–December 2016 were 0.85 (S/Q = 0–1, sample point n = 14,748), 0.38 (S/Q < 0.5, n = 91) and 0.76 (S/Q ≥ 0.5, n = 985).

2.2. Instruments and Measurements

Solar global radiation (Q) and direct radiation (D) were measured at a frequency of 1 Hz using radiometers (model TBQ-4-1 and TBS-2, 322, Institute of Jinzhou, Jinzhou, China) at Qianyanzhou Station, 800 m away from the flux tower [35]. Diffuse radiation (S) is derived from Q-D × cosZ, where Z is the solar zenith angle. PAR was measured using a LI-190SA Quantum Sensor (LI-COR, Inc., Lincoln, NE, USA, relative error < ±5%), except from 1 January to 21 May 2013, when Q and PAR were obtained from the weather station at Qianyanzhou Station. All solar radiation sensors were cleaned every morning and when needed, such as after rain. Detailed explanations of the solar radiation system are reported in [58]. Meteorological parameters (temperature and relative humidity) were measured using a HOBO weather station (Model H21, Onset Company, Bourne, MA, USA) [35]. The measurement frequency was 1 Hz for solar radiation and meteorological parameters.

CO₂ fluxes, together with another set of solar radiation and meteorological parameters, were measured at a 42-m tower starting in 2003 [59]. The eddy covariance technique has been extensively used for the direct measurement of net carbon exchange between the atmosphere and ecosystems [60]. The eddy covariance system included a 3-D sonic anemometer (Model CSAT3, Campbell Scientific Inc., Logan, UT, USA), open-path CO₂/H₂O analyzer (Model LI-7500, Li-cor Inc., Lincoln, NE, USA) and a CR5000 data logger (Campbell Scientific Inc., USA) located 23.6 m above the ground to measure CO₂ flux at the Qianyanzhou subtropical evergreen conifer plantation [59], and observational data (NEE, GPP) measured at 23.6 m from 1 January 2013 to 31 December 2016 were used in this study.

Additionally, BVOC emission fluxes were measured using a relaxed eddy accumulation (REA) system, which was installed on a platform 23 m off the ground in the flux tower [35]. The REA system was composed of a 3-dimensional sonic anemometer (RM Young, Traverse City, MI, USA, Model 81000), a data logger (Campbell Scientific, Logan, UT, USA, Model CR1000) and a data acquisition and control unit and was used to collect air samples in stainless steel cartridges filled with Tenax GR and Carbograph 5TD (Markes International Inc., Gold River, CA, USA). The cartridges were analyzed in the laboratories at the National Center for Atmospheric Research (NCAR) in Boulder, CO, USA and at the Institute of Atmospheric Physics, Chinese Academy of Sciences (IAP, CAS) in Beijing, China. The procedures for sample analysis by gas chromatograph equipped with flame ionization detector (GC-FID) at the Beijing laboratory (IAP, CAS) were similar to those described by Greenberg et al. [32,35]. BVOC emission fluxes were measured from 22 May 2013 to 4 January 2016: air samples of BVOCs were collected around noon or five times with 3 hours interval per day from sunrise to sunset, and a more detailed introduction to BVOC studies is reported by Bai et al. [35]. The measurement frequency was 10 Hz over 30 min for CO₂ flux, NEE and BVOC data collection.

2.3. Flux Data Processing and Data Selection

To ensure reliable processing of flux data, ChinaFLUX developed standard methodologies for processing flux data, including data quality control; coordinate rotation; Webb, Pearman and Leuning density correction (WPL) correction (correction for density fluctuations); canopy storage calculation; gap filling; and flux component partitioning [61]. These methods were applied to process the flux data from the Qianyanzhou site using the software program MATLAB R2014a (MathWorks Inc., Natick, MA, USA). The raw data were corrected using double coordinate rotations [62] and WPL correction [63]. Abnormal or missing data were mainly caused by rainfall, dew, frost, power failure or instrumental malfunction. The Lloyd–Taylor and Michaelis–Menten equations were used for gap filling of carbon flux at nighttime and in daytime, respectively [64,65]. The observed nighttime carbon flux equals ecosystem respiration (Re), which represents the CO₂ released by vegetation and soil respiration. The Lloyd–Taylor equation was also used to extrapolate the daytime Re. The directly observed daytime carbon flux was net ecosystem CO₂ exchange (NEE), whose absolute value equals net ecosystem productivity (NEP, NEP = −NEE). After the calculations of daytime Re, GPP was determined by the equation GPP = Re + NEP and regarded as the observed GPP in this study. Based on the continuous daytime and nighttime datasets, GPP was obtained [66]. The threshold method was used to remove spurious environmental data, and the data gaps were filled using linear interpolation or mean diurnal variation methods [67].

2.4. Empirical Model of Gross Primary Production (EMGPP)

Based on PAR energy balance at the canopy level, an empirical model of gross primary production (EMGPP) considering 3-factor and 2-factor situations was developed. It was a further application of an empirical model simulating BVOC emission fluxes (EMBEs) that was previously developed for this conifer plantation, a temperate forest and a subtropical bamboo forest in China [27]. This EMBE model can reasonably estimate BVOC emissions for the above forests in China as assessed by different evaluations. At the Qianyanzhou conifer plantation, the measured mean and maximum emission fluxes in 2013 were 0.137 and 1.610 mg m⁻²h⁻¹ for isoprene, respectively, and 0.474 and 2.711 mg m⁻²h⁻¹ for monoterpenes, respectively. Monoterpene emission fluxes dominated the total BVOC emission fluxes with a contribution of 71.6% of BVOC emission fluxes, displaying clear diurnal variations (higher in the afternoon and lower in the early morning and late evening) and seasonal variations (higher in summer and lower in winter) [35]. Therefore, considering the assuming similarities between the processes and PAR utilization on BVOC emissions (e.g., isoprene production) [68] and GPP, further application of EMBE was investigated, i.e., similar PAR transfer in the atmosphere (absorption, scattering), and different PAR attenuation by BVOCs or CO₂ in the atmosphere, i.e., only one difference in the description of BVOCs and CO₂ term (or GPP term). More specifically, an empirical model development of GPP is reported as follows.

The PAR energy balance on a horizontal plane above the canopy level is calculated using three terms (the GPP term and photochemical and scattering terms, representing CO₂ attenuation, GLP attenuation (except CO₂) and scattering of PAR, respectively) and described as:

PAR = A₁e^−0.1aGPPtm × cos(Z) + A₂e^−kWm × cos(Z) + A₃e^−S/Q + A₀

(1)

Then an empirical model of GPP was developed following procedures similar to those previously used to describe the BVOC emissions at this conifer plantation as well as other forests in China [27,69,70], i.e., a further application of the empirical model of BVOC emissions (EMBE):

e^−0.1aGPPtm × cos(Z) = B₁PAR + B₂e^−kWm × cos(Z) + B₃e^−S/Q + B₀

(2)

where in the GPP term, e^−0.1aGPPtm, a is the attenuation coefficient for CO₂ in the atmosphere (set to 1 mg CO₂⁻¹ m² s⁻¹), m is the optical thickness of the air mass in the center of the averaging window, GPP is the hourly gross primary production in the sampling period (mg CO₂ m⁻² s⁻¹), t (hour, converted from 0.5 h in CO₂ flux) is the sampling period and 0.1 is a normalizing coefficient for GPP and yields e^−0.1aGPPtm × cos(Z) < 1. According to the Beer–Lambert law, the GPP term was expressed as e^−0.1aGPPtm. In the photochemical term, e^−kWm, W = 0.021E × 60, where E is the mean water vapor pressure (hPa) at ground level during the sampling period and k is the average absorption coefficient of water vapor in the wavelength range of 0.70–2.845 μm. In the scattering term, e^−S/Q, S/Q is a scattering factor that objectively describes the relative quantities of GLPs in the atmosphere, including clouds, haze and rain [58]. Cos(Z) = 1/m and Z is the solar zenith angle (degree). A_i and B_i are coefficients that are determined empirically by using an hourly observational dataset of solar radiation, water vapor pressure and GPP.

To reduce the influence of errors in the carbon flux (i.e., NEE, mg CO₂ m⁻² s⁻¹) and solar radiation measurements in the EMGPP model development, the following data were not used: (1) NEE larger than two times the standard deviation calculated from all individual measurements and the corresponding GPP and (2) global solar radiation with a solar zenith angle greater than 75°.

The EMGPP model was developed utilizing a dataset that includes GPP, with PAR, photochemical and scattering terms (3-factor model, Equation (2)). Considering the practicality of the EMGPP model when direct or diffuse solar radiation is not available, a 2-factor model was developed (i.e., without the scattering term in Equation (2)).

e^−0.1aGPPtm × cos(Z) = C₁PAR + C₂e^−kWm × cos(Z) + C₀

(3)

Therefore, the EMGPP models were developed for these two conditions (3-factor and 2-factor).

3. Results

3.1. Model Development and Evaluation for the Subtropical Coniferous Forest

Because PAR transfers in the atmosphere and the interactions between PAR and atmospheric substances are different under different atmospheric conditions, two atmospheric conditions (expressed as the scattering factor, S/Q) were studied, i.e., S/Q < 0.5 and S/Q ≥ 0.5, representing a relatively clean and sunny atmosphere (low GLP loads, e.g., low quantities of clouds and aerosols) with high solar radiation and air temperature and a misty and cloudy atmosphere (high GLP loads, e.g., large amounts of cloud, high aerosol loads, rain) with low solar radiation and air temperature, respectively.

3.1.1. EMGPP Models Using 3-Factor and 2-Factor for S/Q < 0.5 Conditions

Generally, GPP, NEE and BVOC emissions in this coniferous forest show clear diurnal, seasonal and interannual variations and vary with key driving factors and atmospheric conditions, including PAR, air and soil temperature, precipitation, water vapor pressure, vapor pressure deficit (VPD), soil water content, cloud density, drought and other factors [35,59,71,72,73] and are higher around noon and in summer. Daily emission peaks occur around noon.

Based on the above data criteria in Section 2.4, 91 measurements (n = 91 for S/Q < 0.5) from 22 May 2013 to 31 December 2014 were selected in the model developments for using the 3-factor and 2-factor methods. These valid data were used together with the mean water vapor pressure at ground level (E) and the solar radiation measurements (PAR, S and Q) averaged over the duration of each NEE measurement to determine the coefficients B_i of the 3-factor and C_i of the 2-factor EMGPP models.

Table 1. Coefficients of the EMGPP model (3-factor and 2-factor) (Equation (2)) determined using carbon flux and solar radiation measurements under low GLPs (S/Q < 0.5) in the Qianyanzhou conifer plantation during 2013–2014, statistical metrics, i.e., coefficient of determination (R²), average and maximum of the absolute relative bias (δ_avg, δ_max (%)), normalized mean square error (NMSE), and standard deviations of calculated and observed fluxes (σ_cal and σ_obs), together with mean bias errors (MAD, mg CO₂ m⁻² s⁻¹ and %) and root mean square errors (RMSE, mg CO₂ m⁻² s⁻¹ and %) (n = 91).

Model	B1	B2	B3	B0	R2	δavg	δmax	NMSE	σcal	σobs	MAD		RMSE
											(mg CO2 m−2 s−1)	(%)	(mg CO2 m−2 s−1)	(%)
3 factor	0.029	0.992	−0.087	0.085	0.999	9.96	44.04	0.013	0.14	0.13	0.06	9.46	0.07	11.33
Model	C1	C2		C0	R2	δavg	δmax	NMSE	σcal	σobs	MAD		RMSE
											(mg CO2 m−2 s−1)	(%)	(mg CO2 m−2 s−1)	(%)
2 factor	0.027	1.014		0.025	0.999	10.25	29.49	0.016	0.14	0.13	0.06	10.12	0.08	12.55

shows the results and statistical metrics, including coefficient of determination (R²), average and maximum of the absolute relative bias, δ_avg and δ_max (%), (δ = |y_cal − y_obs|*100/|y_obs|, where y_cal and y_obs are calculated and observed GPP), normalized mean square error (NMSE = $\overline{{(y_{cal}-y_{obs})}^2}/({\overline{y}}_{cal}\times {\overline{y}}_{obs})$) [74], standard deviations of calculated and observed fluxes (σ_cal and σ_obs), mean absolute deviations (MAD, in mgCO₂ m⁻² s⁻¹ and percentage of mean measured value, %), and root mean square errors (RMSE, in mgCO₂ m⁻² s⁻¹ and in percentage of mean measured value).

Figure 1 displays the calculated hourly GPP using 3-factor and 2-factor EMGPP models and measured GPP for S/Q < 0.5 conditions. Figure 2 shows a scatter plot of calculated versus measured hourly GPP. Generally, the 3-factor and 2-factor EMGPP models agree with the measurements for S/Q < 0.5 conditions.

3.1.2. EMGPP Models Using 3-Factor and 2-Factor for S/Q ≥ 0.5 Conditions

It is necessary to develop an EMGPP model for high GLP conditions (i.e., S/Q ≥ 0.5). The data criteria were similar to those for S/Q < 0.5, but the only difference was for S/Q, i.e., S/Q ≥ 0.5, and all other corresponding observational data were used in EMGPP model development for S/Q ≥ 0.5. Similarly, the simulated and observed hourly GPP, along with statistical metrics, is reported in

Table 2. Same as Table 1 but for S/Q ≥ 0.5 (n = 985).

Model	B1	B2	B3	B0	R2	δavg	δmax	NMSE	σcal	σobs	MAD		RMSE
											(mg CO2 m−2 s−1)	(%)	(mg CO2 m−2 s−1)	(%)
3 factor	0.012	1.128	−0.092	0.072	0.994	15.52	47.91	0.035	0.228	0.226	0.12	15.21	0.15	18.74
Model	C1	C2		C0	R2	δavg	δmax	NMSE	σcal	σobs	MAD		RMSE
											(mg CO2 m−2 s−1)	(%)	(mg CO2 m−2 s−1)	(%)
2 factor	0.005	1.171		0.027	0.993	16.20	67.15	0.037	0.235	0.226	0.13	15.66	0.16	19.35

. In general, the calculated hourly GPP was also in line with these measured values, but with larger estimation errors (δ_avg, NMSE, MAD and RMSE) compared to those for S/Q < 0.5, which was caused by a larger observational error (σ_obs) of 0.23 mg CO₂ m⁻² s⁻¹ at S/Q ≥ 0.5 compared to 0.13 mg CO₂ m⁻² s⁻¹ at S/Q < 0.5. Figure 3 shows a scatter plot of calculated versus measured hourly GPP under all sky conditions for S/Q ≥ 0.5.

3.2. Evaluation of the EMGPP Models under All Sky Conditions

3.2.1. Evaluation of EMGPP Models for S/Q < 0.5 Conditions

First, we investigated the performance of 3-factor and 2-factor EMGPP models for GPP estimates from 22 May 2013 to 31 December 2016. The hourly GPP calculated using the coefficients in Table 1 and Table 2 was compared with the observations for S/Q < 0.5. To thoroughly investigate performance in detail, we also show the results for the periods of 2013–2014, 2015, 2016 and 2013–2016. The results in

Table 3. Hourly GPP simulations in 2013–2016 using 3-factor EMGPP for S/Q < 0.5. The statistical metrics are the same as in Table 1, along with the ratios of average calculated versus to the observed GPP during the different time periods (cal/obs).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013–2014	83.34	0.365	0.382	0.247	1.01	873	0.259	42.27	0.510	83.08
2015	62.61	0.283	0.334	0.258	1.03	445	0.265	41.60	0.515	80.90
2016	90.21	0.460	0.389	0.301	1.07	322	0.335	53.29	0.580	92.23
2013–2016	78.41	0.357	0.369	0.260	1.03	1665	0.274	44.08	0.378	60.69

and

Table 4. Same as Table 3 but using the 2-factor EMGPP model (S/Q < 0.5).

Time Period	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
							(mg CO2 m−2 s−1)	(%)	(mg CO2 m−2 s−1)	(%)
2013–2014	76.46	0.333	0.369	0.247	1.04	873	0.252	41.03	0.502	81.84
2015	62.44	0.272	0.329	0.258	1.04	445	0.265	41.59	0.515	80.89
2016	97.44	0.508	0.420	0.301	1.11	322	0.361	57.43	0.602	95.74
2013–2016	76.19	0.349	0.368	0.260	1.05	1665	0.276	44.27	0.378	60.65

and Figure 4 show that the estimated monthly GPP using 3-factor and 2-factor models overestimated GPP with a relative bias of 78.41% and 76.19% during 2013–2016, respectively, and their corresponding ratios of the mean calculated hourly GPP to the observed were 1.03 and 1.05. Hourly GPP calculated using the 2-factor model was generally slightly different from that calculated using the 3-factor model. Overall, the 3-factor modeled GPP was somewhat closer to the observed GPP. It is noted that it is not an accurate validation of the EMBE models for these same data in 2013 and 2014 were used/overlapped for model development and testing, and the others without 2013 and 2014 data (e.g., 2015, 2016 and 2015–2016) were used as validations.

It should be noted that hourly, monthly and annual amounts/averages of GPP were calculated starting from June 2013 for 2013 (i.e., the hourly data during January–May 2013 were not considered), resulting in the amount of hourly sampling points (2013–2014, 2015, 2016) being less than that in 2013–2016 in Table 3.

According to the hourly GPP estimates, daily, monthly and annual averages and amounts of GPP were obtained. Monthly amounts of simulated and observed GPP for S/Q < 0.5 are shown in Figure 4 and

Table 5. Same as Table 3 but for simulations of monthly amounts of GPP (mg CO₂ m⁻²) (S/Q < 0.5).

3-F Time Period	δavg	δmax	GPPcal	GPPobs	cal/obs	2-F Time Period	δavg	δmax	GPPcal	GPPobs	cal/obs
2013	51.08	137.92	40.74	36.32	0.84	2013	37.88	105.36	41.21	36.32	0.97
2014	33.62	87.80	21.23	23.43	0.77	2014	30.09	86.06	22.14	23.43	0.83
2015	26.82	56.51	29.15	28.35	0.95	2015	22.28	54.46	29.48	28.35	0.99
2016	36.10	139.35	19.72	18.40	1.16	2016	43.52	139.93	20.36	18.40	1.24
2013–2014	40.05	137.92	28.42	28.18	0.80	2013–2014	32.96	105.36	29.16	28.18	0.88
2013–2016	35.66	139.35	26.21	25.53	0.93	2013–2016	33.19	139.93	26.82	25.53	1.01

. A scatter plot of calculated versus measured monthly amounts of GPP is displayed in Figure 5. The 3-factor and 2-factor empirical models can well capture the GPP monthly and seasonal variations, i.e., higher values in summer and lower values in winter [72]. Generally, the 3-factor model performs better than the 2-factor model. The 2-factor model can also be used when the scattering factor (e.g., direct or diffuse solar radiation) is not available. The estimated monthly GPP using 3-factor and 2-factor models overestimated GPP with relative biases of 35.66% and 33.19% during 2013–2016, respectively.

Similarly, comparisons of annual amounts of GPP are given in

Table 6. Same as Table 5 but for annual amounts of GPP simulations (mg CO₂ m⁻²) for S/Q < 0.5.

3-F Time Period	δavg	GPPcal	GPPobs	cal/obs	2-F Time Period	δavg	GPPcal	GPPobs	cal/obs
2013	12.17	285.17	254.23	1.12	2013	13.46	288.45	254.23	1.13
2014	9.38	254.74	281.11	0.91	2014	5.50	265.66	281.11	0.95
2015	2.83	291.47	283.46	1.03	2015	4.00	294.79	283.46	1.04
2016	7.20	216.94	202.36	1.07	2016	10.68	223.98	202.36	1.11
2013–2014	0.85	539.91	535.34	1.01	2013–2014	3.51	554.11	535.34	1.04
2013–2016	2.66	1048.32	1021.16	1.03	2013–2016	5.06	1072.87	1021.16	1.05

, and better estimates were also performed, e.g., the ratios of annual amounts of GPP calculated as observed were 1.03 (ranging from 0.91 to 1.12) and 1.05 (0.95–1.13). Additionally, the relative errors were 2.66% (7.20–12.17%) and 5.06% (4.00–13.46%), respectively, for using 3-factor and 2-factor models

3.2.2. Evaluation of the EMGPP Model for S/Q ≥ 0.5

Similar to the above GPP evaluations for S/Q < 0.5, hourly GPP was also computed using the 3-factor and 2-factor models and compared against the observations from 22 May 2013 to 31 December 2016 for S/Q ≥ 0.5. Compared to measured GPP, the modeled hourly GPP was overestimated by 41% (range from 28% to 54%) in 2013–2016 (

Table 7. Same as Table 3, but for hourly GPP estimates in 2013–2016 using 3-factor EMGPP (S/Q ≥ 0.5).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mg CO2 m−2 s−1)	(%)	(mg CO2 m−2 s−1)	(%)
2013–2014	143.80	0.582	0.388	0.293	1.42	5759	0.295	70.83	0.543	130.46
2015	156.93	0.653	0.387	0.287	1.54	3617	0.332	79.89	0.576	138.75
2016	156.01	0.626	0.368	0.312	1.28	3552	0.349	73.65	0.591	124.73
2013–2016	150.73	0.611	0.382	0.298	1.41	13031	0.320	73.93	0.402	92.83

) when using the 3-factor model and 38% (26% to 50%) (

Table 8. Same as Table 7 but for hourly GPP estimates using 2-factor EMGPP for S/Q ≥ 0.5.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mg CO2 m−2 s−1)	(%)	(mg CO2 m−2 s−1)	(%)
2013–2014	137.05	0.522	0.368	0.293	1.39	5759	0.282	67.83	0.531	127.66
2015	148.56	0.557	0.366	0.287	1.50	3617	0.306	73.83	0.554	133.34
2016	146.65	0.576	0.353	0.312	1.26	3552	0.328	69.25	0.573	120.95
2013–2016	143.08	0.545	0.364	0.298	1.38	13,031	0.302	69.71	0.375	86.71

) when using the 2-factor model. In contrast to the GPP estimates at S/Q < 0.5, the calculation errors (δ_avg, σ_cal, MAD, RMSE) at S/Q ≥ 0.5 increased to different extents, which was mainly caused by the larger and increased observational errors in GPP, solar radiation and meteorological variables. For example, the standard errors (σ_obs) were 0.298 and 0.260 mg CO₂ m⁻² for observed GPP at S/Q ≥ 0.5 and S/Q < 0.5 during 2013–2016, respectively, increasing by 15% at S/Q ≥ 0.5 compared to that at S/Q < 0.5. Similarly, the standard error of the observed S/Q, 0.144 at S/Q ≥ 0.5 also increased compared to 0.075 at S/Q < 0.5 during 2013–2016, implying the high atmospheric substance (e.g., aerosols, clouds, fog) loads and the complicated scattering processes in PAR transfer lead to large errors in observations and calculations at S/Q ≥ 0.5 conditions. In addition, the averages of the observed hourly PAR values during 2013–2016 were 3.641 and 1.540 MJ m⁻² for S/Q ≥ 0.5 and S/Q < 0.5, respectively.

The monthly amounts of the calculated and observed GPP for S/Q ≥ 0.5 are displayed in Figure 6 and

Table 9. Same as Table 5 and for simulations of monthly amounts of GPP (mg CO₂ m⁻²) (S/Q ≥ 0.5).

3-F Time Period	δavg	δmax	GPPcal	GPPobs	cal/obs	2-F Time Period	δavg	δmax	GPPcal	GPPobs	cal/obs
2013	43.62	77.66	169.86	134.17	1.19	2013	42.62	74.40	169.71	134.17	1.19
2014	51.87	92.51	183.79	121.52	1.37	2014	50.90	90.06	178.56	121.52	1.31
2015	50.20	95.21	193.13	125.12	1.45	2015	48.70	89.28	187.62	125.12	1.39
2016	36.99	65.80	179.30	140.17	1.20	2016	39.15	70.01	175.86	140.17	1.16
2013–2014	48.83	92.51	178.66	126.18	1.30	2013–2014	47.85	90.06	175.30	126.18	1.26
2013–2016	45.91	95.21	182.87	129.79	1.32	2013–2016	45.66	90.06	178.89	129.79	1.27

. Overall, most monthly GPP estimates agreed with the observed values within 2 times the standard deviations of the observed values. Monthly amount GPP estimates overestimated by 32% (the range of the ratio of calculated to observed GPP was 0.28–1.90) and 27% (0.29–1.90) during 2013–2016 using 3-factor and 2-factor models, respectively. Both monthly amount GPP simulations using 3-factor and 2-factor models can also easily capture seasonal variations (higher in summer and lower in winter) and have similar GPP values, indicating that both models can be used to compute GPP for S/Q ≥ 0.5. A scatter plot of calculated versus measured monthly amounts of GPP is given in Figure S1 (In Supplementary Material).

Generally, the 3-factor model generates a better simulation than the 2-factor model. The 2-factor model can also be selected. Both the 3-factor and 2-factor models overestimated GPP with a relative bias of 45.91% and 45.66% from 2013–2016, respectively.

The comparisons of annual amounts of GPP are given in Table S1 (In Supplementary Material), and EMGPP models generated reasonable simulations. The ratios of annual amounts of GPP calculated versus observed were 1.41 (ranging from 1.26 to 1.54) and 1.38 (1.25–1.50) using 3-factor and 2-factor models, respectively. It is evident that the estimation errors using 3-factor and 2-factor models increased the most compared to those of S/Q < 0.5 conditions, which is understandable and caused by the increased measurement errors accompanying high GLP loads (including aerosols, clouds, fogs, dews and smogs).

3.3. GPP Simulations and Measurements under All Sky Conditions during 2013–2016

To fully investigate GPP levels and variations under all sky conditions, all hourly GPP was calculated by the 2-factor and 3-factor EMGPP models and combined coefficients for S/Q < 0.5 and S/Q ≥ 0.5 (all hourly GPP during 1 January 2013–21 May 2013 was calculated using the 2-factor EMGPP model, as the direct radiation and S/Q values were not available). All statistical metrics are shown in Tables S2 and S3 (in Supplementary Material) for the 3-factor and 2-factor models, respectively. Under all sky conditions, hourly GPP estimates (including S/Q < 0.5 and S/Q ≥ 0.5) were generally larger than the observed values, and simulation errors (δavg, σ_cal, MAD, RMSE) were larger than those for S/Q < 0.5 conditions and slightly smaller than those for S/Q ≥ 0.5 conditions. Specifically, the mean simulated hourly GPP was overestimated compared to observed GPP by 31% (ranging from 15% to 43%) from 2013–2016; RMSE values were 0.398 mg CO₂ m⁻² s⁻¹ and 87.60% when using the 3-factor model, and the corresponding values were 29% (17%-39%), 0.378 mg CO₂ m⁻² s⁻¹, 83.23% when using the 2-factor model. Generally, the EMGPP models showed similar performance over the model development and nonmodel development periods. Thus, the EMGPP model can be applied to simulate GPP under all skies.

Likewise, the mean simulated daily GPP overestimated the observed GPP by 31% during 2013–2016; RMSE values were 3.55 g C m⁻² and 87.60% when using the 3-factor model, and the corresponding values were 29%, 3.43 g C m⁻² and 83.23% when using the 2-factor model.

Monthly amounts of calculated and observed GPP for S/Q = 0–1 are displayed in Figures S2 and S3 (in Supplementary Material) and Table S4 (in Supplementary Material). Most monthly amounts of GPP estimates showed reasonable agreement with the observations within 2 times the standard deviations of the observed. For using 3-factor and 2-factor models: (1) The estimated monthly amounts of GPP were overestimated compared to observed GPP by 22% (the range of the ratio of calculated to observed GPP was 0.28–1.90) and 18% (0.23–1.86) from 2013–2016, respectively. (2) The R² values between the estimated and observed monthly amounts of GPP were 0.86 and 0.87, respectively. (3) Monthly GPP estimates using 3-factor and 2-factor models can also capture seasonal variations and have similar GPP values when S/Q = 0–1, indicating that both models can be applied to calculate GPP under all sky conditions. (4) The annual mean GPP estimates overestimated the observed GPP by 31% (ranging from 15% to 43%) and 29% (17–39%), respectively.

Annual amounts of GPP over different periods in 2013–2016 were calculated and compared with the observed values (Table S5, in Supplementary Material).

The EMGPP models also displayed reasonable performance; the ratios of the annual amounts of GPP calculated versus observed were 1.41 (ranging from 1.15 to 1.43) and 1.29 (1.17–1.39) when using the 3-factor and 2-factor models, respectively, that is, the empirical models still overestimated the annual amounts of GPP by approximately 41% (15–43%) and 29% (17–39%), respectively, while the 2-factor model presented a slightly better performance. Generally, the model performance for S/Q = 0–1 was better than that for S/Q ≥ 0.5 and worse than that for S/Q < 0.5.

Yale Interactive Terrestrial Biosphere (YIB) estimations reproduce measured GPP seasonality (correlation coefficient R > 0.9) with low biases (−21–15%) at 5 out of 10 sites in China. The calculated GPP shows reasonable seasonality in Qianyanzhou Forest, although the YIB model overestimates its magnitude by 71% [75]. The mean RMSE for GPP is approximately 3.21 g C m⁻² day⁻¹ for 57 evergreen coniferous forest (ENF) sites calculated using the YIB model [76], which is in good agreement with the values of 3.55 and 3.44 g C m⁻² day⁻¹ in 2013–2016 obtained when using 3-factor and 2-factor models, respectively (n = 1320). The YIB model yields mean biases of GPP of −33–50% (equivalent to ratios 0.7–1.5) for 117 out of 195 sites [76]. The ratios are smaller than 1.50 (ranging from 1.15–1.43 and 1.17–1.39 for the 3-factor and 2-factor models, respectively) under different conditions in this study (Table S5, in Supplementary Material). It should be mentioned that the YIB model shows good GPP performance in model evaluations among 17 dynamic global vegetation models (DGVMs) [77].

Specifically, YIB simulations show good performance at flux sites in China [75,76]. The predicted GPP shows reasonable seasonality at an EBF site (CN-Din) in the south, although the model overestimates its magnitude by 71%. For the remaining 4 sites, the model underestimates GPP by 49–90%. The lowest correlations (R < 0.6) are found at two northern grassland sites (CN-Sw2 and CN-Du3) because simulations yield near-zero values in winter, while observations still show high values. Additionally, the simulated biases of annual GPP using the original Biome-BGC model and those from data assimilation are −21.76% and −1.62%, respectively [78]. The RMSE values for estimated GPP at 8 days in the Qianyanzhou coniferous forest were 13.41–20.44 g C m⁻² and 35.35–53.88% between the VPM (vegeration photosynthesis model), EC-LUE (eddy covariance-light use efficiency model), TG (temperature and greenness model), and MODIS-PSN (moderate resolution imaging spectroradiometer-photosynthesis) models compared to the observed GPP [79].

3.4. Sensitivity Analysis of GPP

The sensitivity/response of GPP to changes in its influencing factors (i.e., PAR, E and S/Q changed by −80%, −40%, −20%, +20%, +40%, +80%, +160%) was calculated using the 3-factor EMGPP model. Under realistic atmospheric conditions, the change in GPP due to a change in each factor was studied while keeping all other factors at their original levels. In this sensitivity test, the data in developing the EMGPP models were used for the conditions of S/Q < 0.5 and S/Q ≥ 0.5. Their samples were 91 and 985, respectively. The mean changes in GPP (%) in response to the change in each factor (%) are shown in Figure S4 (in Supplementary Material) and Table S6 (in Supplementary Material) (taking the change rates of 20% in each influencing factor as an example).

It is evident that GPP at this conifer plantation was more sensitive to PAR than to water vapor and S/Q and more sensitive to water vapor than to S/Q at low S/Q, corresponding to better/cleaner atmospheric conditions (high PAR and air temperature and low GLP loads, mainly aerosols, clouds and smog). At high S/Q corresponding to the worse/smoggy atmospheric conditions (low PAR and air temperature, and large GLP loads), GPP was most sensitive to E (meaning that the water supply/drought is a very important and predominant factor in controlling/influencing GPP under most atmospheric conditions) than to PAR and S/Q and was slightly more sensitive to PAR than to S/Q. This result reveals that GPP would decrease significantly with only the increase in PAR at low S/Q (e.g., under the conditions of stomatal closure at noon in the summer) and a much smaller GPP decrease at high S/Q (Table S5, in Supplementary Material). The increase or decrease in GPP with the decrease or increase in PAR reveals a mechanism by which the PAR energy is absorbed and utilized by plants, while all other factors remain at their original levels. GPP will increase similarly to the increase in water vapor under both S/Q conditions, meaning that sufficient water supplies are beneficial to plant growth and ecosystem carbon uptake. GPP will increase with a decrease in S/Q, implying that the decrease in scattering GLPs in the atmosphere results in an increase in diffuse PAR arriving at the canopy, followed by an increase in GPP, which corresponds to enhanced diffuse radiation promoting carbon assimilation [80,81,82,83,84,85,86], and increased cloudiness reduces total radiation, which decreases GPP [87]. The mechanism of and relationship between scattered PAR (or scattered GLPs, S/Q) and GPP are also manifested when higher GLP loads (i.e., S/Q) cause a larger decrease in GPP (compared to lower GLP loads) when S/Q changes at the same rate, and vice versa (Table S5, in Supplementary Material). The responses of GPP to its influencing factors strongly depend on the levels of influencing factors, e.g., larger responses of GPP to PAR at high PAR (4.57 μmol m⁻² s⁻¹, S/Q < 0.5) and to S/Q for high S/Q (0.76, S/Q ≥ 0.5), along with closer responses of GPP to E at close E (~28.0 hPa) under S/Q < 0.5 and S/Q ≥ 0.5 conditions. The average S/Q value of 0.85 under all sky conditions in 2013–2016 indicated that the atmosphere in the Qianyanzhou region was dominated by high GLP loads, i.e., S/Q ≥ 0.5 (Section 2.1). Under these atmospheric conditions (common situations), (1) GPP will increase with climate wetting and decrease with climate drought, and climate drought will lead to a larger influence on GPP increase (i.e., enhancement in carbon fixation capacity) than on GPP loss caused by climate wetting; (2) GPP would increase with the cleaner atmosphere (decrease in S/Q), meaning that the carbon fixation capacity of the forest would benefit from the cleaner air and improved air quality (low GLP concentrations). It should be noted that the response of GPP to the decrease of each factor was slightly larger than to the increase of each factor. Therefore, a cleaner atmosphere and better air quality will be favorable for more carbon fixation in plants, i.e., carbon emission peaks and carbon neutrality will lead to more GPP and slow climate warming. Furthermore, air pollution control and carbon peak and carbon neutrality have beneficial and positive effects on nature and human beings. Similar results to cobenefits of carbon neutrality and air pollution reduction have been reported [88].

In addition, the sensitivity results can also be used as references to evaluate the uncertainties of the empirical model, including with the changes in each factor or their combinations.

4. Discussion

4.1. EMGPP/CO₂, PAR and Related Mechanisms

The EMGPP model expresses PAR transfer and usage in atmosphere–vegetation–land associated with the PAR attenuation of the equivalent CO₂ for producing GPP, absorption and usage of PAR due to GLPs (except CO₂), and GLPs scattering.

In brief, PAR absorption and indirect utilization by GLPs through chemical and photochemical reactions mainly include NO₂ + hv (λ > 420 nm) → NO₂*, NO₂* + H₂O → HONO + OH [89], and BVOCs + OH + O₃. The photochemical term represents the total PAR utilization by all atmospheric GLPs (except CO₂) associated with NO₂* and OH radicals. In Equations (1) and (2), the photochemical term, the scattering term and the complex interactions between PAR and the three terms (absorbing, scattering and BVOC terms, which are replaced by the GPP term in this study) are fully explained in previous papers [27,40,70]. In more specific, the chemical compositions (O₃, NOx, SO₂, HCHO, etc.) react with BVOCs and anthropogenic VOCs (AVOCs) through OH radicals and H₂O; their roles are implicitly described in the photochemical term in EMBE and EMGPP. If the individual roles of these species in GPP and BVOCs need to be studied quantitatively, the specific constituents (e.g., O₃) can be added and explicitly described in EMGPP and EMBE, following the previous successful method [40,90]. The fact that high O₃ results in GPP decline in general [91,92] is an important issue, and the authors suggest studying it by adding an O₃ term to EMGPP in the future.

The EMBE model can reasonably simulate BVOC emission fluxes over different time scales (diurnal to annual) and in different situations (e.g., single measurement and hourly average in the measurement campaign) [27]. Although the BVOC concentrations (usually in ppb or ppt) were very low in the atmosphere, their roles in PAR and interactions with other factors (e.g., photochemical and scattering) can be readily captured and represented in photochemical and scattering terms. Similarly, CO₂ concentration (several hundred ppm) is much higher than BVOCs in the atmosphere and its role can also be captured by using the PAR balance method and expressed in their related terms.

PAR balance and EMGPP (Equations (1) and (2)) express the multiple interactions of PAR, CO₂ and other GLPs associated with the main energy processes of CO₂, absorbing and scattering GLPs, i.e., solar energy-atmospheric substances. Thus, the objective energy relationships were quantified and then applied to compute GPP.

The signs of coefficients A_i are −, +, +, − in Equation (1) (corresponding to B_i +, +, −, + in Equation (2)), demonstrating that GPP resulting from the presence of PAR is attenuated by the equivalent CO₂ in the atmosphere. The negative sign of A₁ reveals that this part of the energy is captured from the photochemical term/system through the interactions between CO₂ molecules and PAR photons, total PAR direct absorption and indirect consumption by all atmospheric GLPs (excluding CO₂). This is similar to the roles of BVOCs in this coniferous forest [27].

4.2. Performance of EMGPP

Data quality plays a more vital role than data quantity in empirical model development and evaluation. For example, though the small observational datasets 8 and 18 were used in EMBE models for isoprene and monoterpenes, respectively [27], the EMBE showed reasonable performance. Compare to EMBE development, more datasets were used in EMGPP development and the GPP simulations (~4 years) were acceptable by similar uncertainties with the observations (e.g., 20–50%) [26] and widely used models (e.g., RMSE).

It is difficult to find hourly results of other models for more comparisons and evaluations of the EMGPP, which may be related to large uncertainties in GPP hourly simulations/processes. This study focuses on the simulation of GPP as a first step, the simulations of respiration and NEE (or the partitioning of GPP), as well as long-term variations in GPP, respiration and NEE since 2003 will be investigated in other studies.

Accurate simulations of PAR and global solar radiation at high GLP loads (S/Q > 0.5) are still significant challenges [93] for GPP calculations, compared to GPP estimates at high with low S/Q levels. More studies in improving estimations of PAR and global solar radiation as well GPP in high S/Q situations are needed.

It should be mentioned that once CO₂ is emitted into the atmosphere, it takes part in the attenuation, absorbing and scattering processes described/captured by GPP (corresponding to total CO₂), photochemical and scattering terms (temperature, water vapor, S/Q factors), respectively, i.e., CO₂ interacts with other atmospheric substances and PAR that is attenuated by the absorption and scattering of GLPs. That means that CO₂/GPP is not controlled by PAR only. This is another reason to study GPP in this study and respiration and NEE in other studies.

Dong et al. [26] reported a relative uncertainty of 20–50% for CO₂ flux measurements by eddy covariance, which depends on high or low CO₂ flux magnitudes. Annual amounts of net ecosystem CO₂ exchange with an error bound of 600 ± 130 g C m⁻² year⁻¹ and 300 ± 180 g C m⁻² year⁻¹ in different forests [24]. The uncertainty in the annual estimation of NEE is 12–32%, as reported by Loescher et al. [25]. Comparisons between annual estimates using traditional ecological models and eddy covariance measurements are in agreement with 5%, 30% and 100% (for NEP) [22,23,24]. The RMSE values of the Breathing Earth System Simulator (BESS) and MODIS 8-day composite GPP products evaluated against flux tower observations are 2.58 (ranging from 1.37–3.86) g C m⁻² d⁻¹ and 2.86 (1.56–4.80) g C m⁻² d⁻¹, respectively [94]. All of the above results are beneficial and useful for understanding the uncertainties in measurements and model simulations.

According to the calculations and evaluations of hourly, daily, monthly and annual GPP and the similar RMSE values for GPP simulations with the widely used models [76,77,78,79], as well as the above uncertainties in CO₂ flux observations and model estimates, EMGPP can be used to simulate GPP in this coniferous forest. As a first step, the application of the EMBE model was suitable for GPP simulations at this site. Based on the previously successful development of EMBE in this forest to other typical forests in China [27], further development of EMGPP for other ecological systems in China, such as the grass ecosystem, is necessary and beneficial, [95,96].

The EMGPP model performances for daily averages during 2013–2016 were similar to those for hourly averages and are not shown for the purpose of saving space. To thoroughly understand the absorbing and scattering processes, it is suggested to measure direct and/or diffuse solar radiation at sites in China and around the world. The 3-factor model provides a more thorough understanding of the processes than the 2-factor model.

There were larger drops in the measured monthly amounts of GPP in July and August 2016, which were found and caused by some problems in CO₂ flux measurements. The simulated hourly and monthly GPP can be used to correct these observed GPP (Figure S2).

The previous model studies are too complicated, and many CO₂ processes/key parameters are simulated along with assumptions, most of which are interconnected and difficult and time-consuming to express independently, whereas, the EMGPP model can well study multiple interactions between solar radiation, atmosphere and vegetation, using a few parameters which are obtained easily from weather and solar radiation stations, requiring fewer assumptions and saving computing time.

4.3. GPP and Its Driving Factors under Different Sky Conditions

Water vapor pressure is an objective factor to represent the dynamic exchange and balance (or vapor pressure deficit, VPD) of water and water vapor among the atmosphere, plants and soil, as well as PAR use by GLPs, i.e., water vapor plays vital roles in not only water circulation but also PAR energy use/transfer (through OH radicals) in atmosphere–plant–land interaction. It should be mentioned that the responses of GPP to the changes in air temperature and relative humidity are extrapolated from and similar to (i.e., positive or negative) those in water vapor, as there is an inherent relation between water vapor and air temperature and relative humidity (T and RH) [97]. Water vapor pressure acts not only as a surrogate for T and RH, but in more significant and objective roles than T and RH in the processes in the atmosphere, vegetation and land. The dominant roles of temperature in GPP or BVOCs can be represented by/reflected in the water vapor factor (such as in the sensitivity test) to some extent. In more detail, water supply and drought play different roles under different sky conditions (e.g., high or low S/Q). For example, precipitation can accelerate GPP in this subtropical Pinus forest [72,98]. Further investigations about water supply (precipitation, drought, etc.) are still needed.

At the canopy interface, chloroplasts absorb the incoming PAR to convert CO₂ and H₂O and decompose water into [H] and O₂. The attenuated PAR by CO₂ (equivalent to GPP and described by GPP term) in the atmosphere was utilized by plants to produce GPP, which is described by the GPP term. The GPP increase/decrease with the decrease/increase in PAR clearly indicates this intrinsic relationship of PAR transfer, use and conversion between the atmosphere and plants (sensitivity analysis). GPP increases with the decrease in S/Q (corresponding to the increase in scattered PAR, i.e., the decrease in S/Q leads to the increase in scattering e^−S/Q), and GPP increases with the increase in water vapor/water supply. These results were in good agreement with several model studies and observations, e.g., the increased diffuse radiation results in the enhancement of GPP [80,81,83,84,85,86,87,99,100], indicating that similar mechanisms exist in the ecosystem in different regions and that the EMGPP empirical model can reasonably describe some mechanisms in GPP. The PAR energy point of view deserves further investigation to study complicated GPP processes and multiple interactions in PAR–atmosphere–vegetation–land. In more detail, there are dynamic exchanges of chemical components (e.g., CO₂, BVOCs, water, particles) in the atmosphere–vegetation (by stoma)–soil. It should be also emphasized that there are multiple interactions in GPP and its influencing factors (PAR, E, S/Q) as well as their related processes (absorption and scattering in the atmosphere), not one-direction interactions acting as assumed in the sensitivity test. It is important and beneficial to develop an energy method to reveal the inherent mechanisms/interactions associated with GPP and capture/quantify its main processes. GPP and all its associated factors/processes change synchronously, and PAR-CO₂-other GLPs-biosphere-land is suggested to be studied as a whole system.

4.4. The Relationship between BVOCs and NEE

The contributions of BVOC emissions (including isoprene and monoterpenes) to NEE and GPP (BVOCs/NEE, BVOCs/GPP) were 3.5% and 0.8%, respectively, from 2013–2016, according to the hourly average estimates of BVOC emissions [27] and GPP and the observed NEE. The ratio of BVOCs/NEE varies from a few percent in most forests to 25% in a mopane savanna woodland, Botswana (19°54 S, 23°33 E) [31,32,33,35]. The ratio of BVOCs/GPP is 0.6% at a ponderosa pine plantation [34] and close to but slightly smaller than that at the Qianyanzhou conifer plantation. The ratios of BVOCs/NEE were 3.5% from 2013–2016, which was an updated value compared to the previous 1.5% using the annual NEE in 2003 [27], 4.0% in a ponderosa pine plantation (38°53′ N, 120°37′ W) [34], less than 1.0% in the Changbai Mountains temperate forest (42°24′ N, 128°6′ E) from June-September (except for 1.4% in September) [70] and 0.63% at a poplar short-rotation plantation (51°6′44″ N, 3°51′2″ E), Belgium [101]. In the Inner Mongolia grassland (43°26′–44°08′ N, 116°04′–117°05′ E), isoprene/NEE was 2.8–4.3% during the growing seasons of 2002 and 2003 [69]. It is speculated that BVOCs/NEE is higher in the tropical zone, followed by a decline in subtropical and then temperate zones. However, grasslands may have a larger ratio than forests, and further confirmation is needed. It should be noted that BVOCs are necessary compositions of carbon emitted by vegetation into the atmosphere; although they make up small fractions in forests and grasslands, they play significant roles in O₃ and aerosol formation as well as cloud formation through chemical and photochemical reactions, influencing solar radiation transfer and distribution in the atmosphere and global warming [5,31,38,43,93,102,103]. In addition, as BVOC emissions are significant to carbon balance on a global scale [31,37], total amounts of BVOCs emitted into the atmosphere should be considered in regional and global carbon balance [34].

Strong linear and nonlinear relationships were found between BVOCs/NEE and climate zones for the ecosystems except the Inner Mongolia grassland and the savanna woodland (Figure S5).

The responses of BVOCs (using the REA technique and taking isoprene as a representative) to changes in each driving factor (i.e., PAR, E, S/Q) were generally much larger than those for GPP, which was observed and confirmed in a subarctic wet heath by [104] that warming increases gross ecosystem production, but these increases are much lower than those for BVOCs, and the responses to PAR and E (E also represents air temperature) were much higher than those to S/Q for BVOCs and GPP [27] (Tables S6 and S7, in Supplementary Material). The nonlinear relationship is more reasonable for revealing intrinsic interactions between BVOCs, GPP and their controlling factors described in EMBE and EMGPP models. These mechanisms are the main reasons for the strong relationship between BVOCs/NEE and climate zones, and BVOCs/NEE values are higher in tropical climates with high PAR and water vapor and lower in temperate climates with low PAR and E. Based on this mechanism (BVOCs/NEE is mainly dependent on PAR and water vapor), it is speculated that BVOCs/NEE is the largest in the tropical climate and the smallest in the Arctic/subarctic, and more evidence is needed. A similar mechanism but for a positive and nonlinear relationship between isoprene/CO₂ and leaf temperature ranging from 20 °C to 40 °C is observed [31].

4.5. Interactions between GLPs (CO₂, BVOCs), Solar Radiation and Climate Change and Their Potential Effects

It was found that the more GLPs there are in the atmosphere, the higher the air temperature is, and a strong positive relationship exists in air temperature and scattering factor at three sites (in the Arctic, Antarctic and mid-latitude in the Northern Hemisphere) in long-term annual averages [5]. In the context of multiple interactions of GLPs (GHGs and non-GHGs, air pollutants, etc.) and solar radiation (global, PAR, UV, etc.), an increasing amount of carbon will be fixed in vegetation in a cleaner atmosphere (sensitivity section), i.e., carbon emission peaks and carbon neutrality and better air quality will have positive effects in controlling climate warming and promote the achievement of UN sustainable development goals (SDGs).

Atmospheric GLPs, coming from primary emissions (e.g., NO₂, NO, SO₂, AVOCs, BVOCs and particulate matter) and secondary production (e.g., O₃, secondary organic compounds, black carbon and HCHO) through chemical and photochemical reactions triggered by solar UV and visible radiation, change with solar radiation over horizontal and vertical scales. The interactions in solar radiation and all kinds of absorbing and scattering GLPs are region dependent. Some GLPs absorb solar radiation in different wavelength bands (UV, visible and near-infrared) [38,40,43,93], and others indirectly use solar radiation energy by reacting with these absorbing GLPs, OH radicals, water, water vapor and reactive BVOCs, including C₂H₂O₂, CH₃CO radical, NO₃ radical, OClO, CHOCHO, biacetyl, butenedial, NOCl, and thousands of AVOCs and BVOCs [40,46,89,103,105,106,107]. In addition, pyruvic acid and all three α-dicarbonyl compounds have absorptions in the visible and UV regions [108]. Thus, GLPs and solar radiation interact in 3 dimensions regionally and globally and control climate and climate change. CO₂, GHGs and non-GHGs, along with other GLPs, should be considered together in mitigating global warming, as they are undergoing different absorbing and scattering processes and later contributing to climate change [5,93]. In addition, astronomical factors should also be considered [109].

According to all the above results, i.e., GPP simulations, sensitivity analysis and its revealed mechanisms and BVOCs/NEE, the PAR energy method is a good option and tool to explore and understand the key processes as well as complex interactions/mechanisms between radiation–atmosphere–biosphere–land and to avoid considering a large number of specific processes and their complicated expressions. The application of EMGPP is time-saving, makes fewer assumptions, and uses inputs which are easily obtained at stations. They are the innovations of using energy balance principle.

5. Conclusions

With further application of a previously developed empirical model of BVOC emissions at a subtropical conifer plantation, China, an empirical model of GPP in this plantation was developed for a better understanding of the processes of GPP and the relationships in GPP and its driving factors. In model development using 2013–2014 observations, 3-factor and 2-factor EMGPP models generally performed good simulations for hourly, daily, monthly and annual GPP values. For example, when using a 3-factor model, the hourly GPP overestimated the observed GPP by 9.96% and 15.52%, respectively, for S/Q < 0.5 and S/Q ≥ 0.5 situations, and their corresponding RMSE values were 0.07 mg CO₂ m⁻² s⁻¹ and 11.33% and 0.15 mg CO₂ m⁻² s⁻¹ and 18.74%. Evaluations were conducted for hourly, daily, monthly and annual GPP against observations in 2013–2016 under all skies (i.e., all S/Q conditions), and reasonable agreement was found between simulated and observed GPP. Generally, the modeled GPP values were still larger than the measured values, and their standard deviations were relatively close. The GPP in 2013–2016 under all sky conditions was computed using 3-factor and 2-factor EMGPP models, and the simulated annual amount of GPP overestimated these observed values by 41% and 29%, respectively. The GPP empirical model captured some characteristics, such as diurnal and seasonal variation patterns. The 2-factor empirical model indicated similar performance as the 3-factor model and can be used when diffuse or direct solar radiation is not available.

Sensitivity analysis demonstrated that the responses of GPP to its driving factors (PAR, water vapor and S/Q) evidently changed at different atmospheric GLP loads (i.e., S/Q < 0.5 and S/Q ≥ 0.5, the latter condition being a predominant state of the atmosphere), namely, GPP was more sensitive to changes in PAR than that in E and S/Q at low S/Q, while GPP was most sensitive to the changes in E than to PAR and S/Q at high S/Q. Therefore, atmospheric column GLP contents should be taken into consideration in the processes of GPP and the interactions in PAR–CO₂–other GLPs. Some inherent relations between GPP and scattered PAR, water vapor, and atmospheric scattered substances were investigated extensively, which were confirmed by other observations and model studies.

The contributions of BVOC emissions to NEE (BVOCs/NEE) were small (3.5%) in this coniferous forest and other forests in the world and exhibited a clear latitudinal variation. It is assumed that BVOCs/NEE was the largest in the tropical region and lowest in the Arctic/subarctic region, which was partially confirmed by laboratory measurements and sensitivity analyses. BVOCs/GPP was 0.8% in this coniferous forest. BVOCs, due to their very high reactivity, produce large quantities of new atmospheric GLPs, and they are also a necessary part of carbon in the atmosphere and vegetation; thus, they should be considered in regional and global carbon balance, as well as climate change and global warming. Reducing both carbon emissions and air pollutants has beneficial effects on nature.