Empirical Models of Respiration and Net Ecosystem Productivity and Their Applications in a Subtropical Coniferous Plantation in China

Abstract

Net ecosystem exchange (NEE), solar radiation (including photosynthetically active radiation PAR), and meteorological parameters were measured in a subtropical coniferous plantation in China during 2013–2016. Applying the PAR balance principle at a canopy level and analyzing the observation data, an empirical model of respiration (Re, EMRe) considering 3-factor and 2-factor situations was developed and tested for all sky conditions. Generally, the respiration simulations were in reasonable agreement with the observations for the hourly, monthly, and annual sums of respiration. For example, using 3-factor and 2-factor models, the estimated annual sums of daytime and nighttime respiration in 2013–2016 overestimated that which was observed by about 31% and 26%, respectively. Further applications of EMRe and an empirical model of gross primary production (GPP, EMGPP) developed previously at this site, and an empirical model of net ecosystem productivity (NEP, EMNEP) using 3-factor and 2-factor models were obtained (NEP = GPP-Re) and evaluated for all sky conditions. Generally, the simulations of the hourly, monthly, and annual sums of NEP showed reasonable performances. The estimated NEP values overestimated the observations by 22% and 27% for the hourly sums in 2013–2016 when using the 3-factor and 2-factor models, respectively, and 7% and 12% for annual sums in 2013–2015 (2016 data were not used as the CO₂ flux measurements had some problems in the 2016 summer). The NEP estimations were evidently improved when more factors (e.g., dark respiration) influencing Re were considered in the daytime respiration compared to those without considering these factors. To simplify the numerous and complicated CO₂ processes in the simulations of Re and NEP, the PAR energy method was applied to capture and describe its main processes and energy interactions. The PAR energy method was suitable for studying the energy relationships associated with CO₂ processes and developing empirical models for the simulations of GPP, Re, and NEP. These models were useful tools to investigate the multiple interactions and mechanisms between CO₂, other atmospheric compositions, and PAR. Thus, the energy method is suggested to be applied to carbon balance.

1. Introduction

Global warming (especially in three Polar regions and other representative regions in the 20th century) is reported by the Intergovernmental Panel on Climate Change (IPCC) and several studies have aroused the great concern of scientists and government officials [1,2,3,4,5,6,7]. The global CO₂ concentration has been steadily rising from 1959 to 2021 (https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_annmean_mlo.txt, accessed on 3 August 2023). CO₂ is accounting for 70% of greenhouse gases (GHG) emissions with an atmospheric lifetime of 200 years [8]; therefore, the emissions and increases in GHGs, (mainly CO₂) [9,10] should be reduced in the context of global carbon neutrality. The terrestrial biosphere sequesters 20–30% of global anthropogenic CO₂ emissions [11]; thus, it is an important task to accurately evaluate the regional and global ecosystem carbon balance, to better contribute to achieving the international double carbon goal.

As an important basis, measuring net carbon exchange by forests is essential to obtain reliable carbon balance (e.g., source, sink, exchange) and validate/evaluate the model performances. As a second step, all kinds of models can be developed and then applied to 1) study processes and mechanisms of CO₂ flux, net ecosystem exchange (NEE), net ecosystem productivity (NEP) (NEE= -NEP), respiration (Re) and gross primary production (GPP), 2) simulate the variations of NEP, Re and GPP, 3) understand the responses of NEE, Re, GPP to their driving factors, as well as to climate change [12,13,14,15]. For example, these models include empirical models, sophisticated models [14,16,17,18], and atmospheric inversions [19,20]. The popular empirical models consider several processes (mostly in one direction) that are of limited in understanding the mechanisms. The complicated models can elucidate the specific processes in vegetation, soil, and atmosphere and obtain a better understanding of the processes and their mechanisms. The widely used process-based dynamic global vegetation models (DGVMs) are useful tools to calculate the natural land sink, and the DGVM multi-model mean estimate of the global land sink is in agreement with the global carbon budget residual land sink. Currently, novel and innovative machine learning algorithms offer powerful tools to estimate carbon balance (GPP, NEE, Re, etc.), especially extreme gradient boosting (XGBoost), artificial neural networks (ANN), and GIS data-based studies. These artificial intelligence might meet the challenges in carbon balance. Until now, there are still large differences among models [14,15,18,21,22,23,24,25,26,27,28,29]. Additionally, many hypotheses and assumed parameters were used in these complex models, which were limited by our current understanding of the carbon balance processes [21]. The Eddy covariance (EC) technique is reliable for obtaining CO₂ exchange between the atmosphere and ecosystem, and EC flux data are commonly used for validations [22,23]. There are also large uncertainties in CO₂ flux observations and model predictions [22,23,24,25,26,27,28,29]. On the basis of current progress and urgent needs, the internal processes and external drivers, as well as their relationships between solar energy, atmosphere, vegetation, and land, should be fully studied.

Respiration is an important and integral part in regional and global carbon cycling [30,31]. O₂ uptake and the subsequent metabolism of the reaction product in the photorespiratory process is the second highest metabolite flux in plants behind photosynthesis [32,33]. New information about respiration and GPP to ensure a better partitioning between the components of the CO₂ atmospheric budget is necessary [34]. NEE or NEP is a vital indicator to evaluate ecosystem carbon exchange, whereas the accurate simulation of respiration is still a big challenge [31], as it is very small compared to NEE and GPP, and many factors/processes control daytime and nighttime respirations through different mechanisms.

Three metabolic processes control carbon fluxes in vegetation, namely, photosynthesis, photorespiration, and respiration, which are significant and influenced by many factors (temperature, light, water, relative humidity, O₂ and CO₂ concentrations, etc.). Many photosynthesis models (e.g., photosynthesis–irradiance curves, non-rectangular hyperbolic, or rectangular hyperbolic models) were developed and commonly used in their simulations, but there are few photorespiration models [35,36]. In addition, the Lloyd–Taylor equation was widely used to extrapolate the daytime respiration from nighttime respiration [37,38].

China has diverse ecosystems covering several climate zones. Chinese scientists have implemented many studies on China’s land carbon budgets and made much progress. However, there is still the need to accurately evaluate the terrestrial carbon budgets and reduce the uncertainties, considering the large differences in the magnitude of the land carbon budget among numerous model simulations [15,39,40,41,42,43,44].

In order to simplify the sophisticated description in the respiration and NEP processes and to avoid using too many hypotheses and assumed variables, we further applied the photosynthetically active radiation (PAR) balance principle in the respiration and NEP simulations, based on reasonable simulations of GPP using an empirical model of GPP (EMGPP) and PAR balance method [45].

The objective of this study is to (1) develop an empirical model of respiration (EMRe), and an empirical model of NEP (EMNEP) combined using EMRe and EMGPP developed earlier for a thorough understanding and estimate carbon balance in a subtropical coniferous forest in China; (2) evaluate these two empirical model performances against the observations; (3) simulate the Re and NEP in 2013–2016 under all sky conditions; (4) test the feasibility of the PAR balance method and understand the carbon exchange thoroughly by applying the PAR balance method to simulate the Re and NEP from a previous simulation of GPP.

2. Instrumentation and Methods

2.1. Site Description

The measurements of CO₂ fluxes, solar radiation, and meteorological parameters were carried out from a 45 m tower installed by the Chinese Academy of Sciences (CAS) at the Qianyanzhou subtropical coniferous plantation, Taihe county, Jiangxi province (26°44′48″ N, 115°04′13″ E, 110.8 m) since 2003 [46]. BVOC emission fluxes, solar radiation, and meteorological variables were also measured at this tower from 22 May 2013 to 4 January 2016 [47]. Main trees are pinus massoniana, pinus elliottii, cunninghamia lanceolata, and broadleaf trees, with a coverage of 42.0%, 49.6%, 5.8%, and 2.6%, respectively. The main shrubs are Loropetalum Chinense, Adinandra millettii, and Lyonia Compta. The forest surrounding the tower covers about 90% in 1 km² and 70% in 100 km², and the average canopy height is 18 m. The mean slope of this study region is 2.8–13.5°. Annual precipitation is 1485.1 mm, and annual temperature is 17.9 °C [48]. The annual sum of global solar radiation was 4579 and 4496 MJ m⁻² in 2013 and 2014, respectively, and the annual sum of PAR was 7998 and 7656 mol m⁻² in 2013 and 2014, respectively [47].

2.2. Instruments and Measurements

The EC system including 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 CR5000 datalogger (Campbell Scientific Inc.) was located at 23 and 39 m above the ground to measure CO₂ flux at Qianyanzhou subtropical coniferous plantation [46], and the observational data (NEE, Re and GPP) at 23 m during 1 January 2013 to 31 December 2016 were used in this study.

Solar global radiation (Q), direct radiation (D), and PAR were measured using radiometers (model TBQ-4-1 and TBS-2, 322 Institute of Jinzhou, China; LI-190SA Quantum Sensor, LI-COR, Inc., Lincoln, NE, USA,) and a solar radiation system with a frequency of 1 Hz at the Qianyanzhou Experimental Station (CAS), 800 m away from the flux tower [47]. The diffuse radiation (S) is derived from Q-D×cosZ, Z is the solar zenith angle. Q, D, and PAR were measured routinely, except from 1 January to 21 May 2013, when Q and PAR were obtained from the weather station at the Qianyanzhou Experimental Station. All solar radiation sensors were cleaned every morning and when they are needed. Meteorological parameters, temperature (°C), relative humidity (RH, %), and water vapor pressure (E, hPa) were measured using a HOBO weather station (Model H21, Onset Company, Bourne, MA, USA) [47].

2.3. Flux Data Processing and Usage

ChinaFLUX developed some methods for evaluating the performance of observation system and flux data quality control (e.g., coordinate rotation, WPL correction, canopy storage calculation, nighttime flux correction, gap filling, and flux partitioning [49]), which were used to deal with the flux data in this study. The raw data were processed using double coordinate rotations, Webb, Pearman, and Leuning density correction [50,51]. Two equations were used for gap filling of CO₂ flux in daytime and nighttime [37,52]. The observed daytime carbon flux was NEP. After the calculations of daytime respiration, the observed GPP was calculated by GPP = NEP + Re and the routine daytime and nighttime observations [38]. More detailed introduction (threshold method, data gap filling, etc.) is reported in studies [45,53,54].

On the EMRe model development, the following data were not used so as to reduce the influence of measurement errors of carbon flux (mgCO₂ m⁻² s⁻¹) and solar radiation: (1) NEE larger than two times the standard deviation computed from all individual observations along with these corresponding respirations and GPP values, (2) global solar radiation with Z > 75°.

2.4. Empirical Model of Respiration (EMRe)

Based on PAR balance at canopy level, an empirical model of BVOC emissions (EMBE) was developed for some representative forests (a subtropical coniferous forest, a subtropical bamboo forest, and a temperate forest) in China [55]. Further application of PAR balance principle and EMBE, an empirical model of GPP (EMGPP) was developed in Qianyanzhou subtropical forest with reasonable performance [45], considering the similarities of PAR use in the processes of BVOC emissions and GPP, in which CO₂ (GPP) was replaced by BVOCs in the EMGPP. Furthermore, these experiences were used in empirical model developments of respiration and NEP in this study.

Briefly, EMGPP (Equation (1)) considers the following three main processes: PAR attenuation by CO₂ in the atmosphere (GPP term); PAR attenuation by atmospheric gases, liquids, and particles (GLPs, except CO₂, photochemical term); PAR scattering by GLPs (scattering term [45,55]). At a canopy level, these three processes contribute to PAR balance. After adjustment, the EMGPP is described as follows:

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

(1)

where a = 1 mgCO₂⁻¹ m² s⁻¹ is an attenuation coefficient for CO₂ in the atmosphere, m is the optical thickness, GPP is the hourly GPP in the sampling period (mgCO₂ m⁻² s⁻¹), sampling period t = 1 h, and 0.1 is a normalizing coefficient. W = 0.021E × 60, E is the average water vapor pressure (hPa) at the ground in the sampling time period, and k is the mean absorption coefficient of water vapor in 0.70–2.845 μm. Scattering factor (S/Q) represents the atmospheric GLP amounts [47]. Z is the solar zenith angle (degree), and Cos(Z) = 1/m. Coefficients A_i are determined using observational dataset (solar radiation, water vapor pressure, GPP, etc.).

It is known that GPP = Re + NEP, and Re is ecosystem respiration. Then, it is necessary to calculate ecosystem respiration. Firstly, an empirical model of daytime respiration (mgCO₂ m⁻² s⁻¹, Equation (2), 3-factor) (under conditions of light) was developed and named as Re (LR) method.

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

(2)

where, we set b = 1, m = 1 and cos(Z) = 1 for daytime estimation. Coefficients B_i (i = 1, 2, 3) and B₀ were determined by analyzing observational data (data criteria, Section 2.3) in the daytime (Re, solar radiation, water vapor). e^−0.1bRetm × cos(Z) > 1 was used in Re model development. The 3-factor and 2-factor EMRe models were developed for two atmospheric situations (S/Q < 0.5 and S/Q ≥ 0.5, representing low atmospheric GLPs, high solar radiation and air temperature, and high GLPs, low solar radiation and air temperature, respectively), which were similar to EMGPP and EMBE models [45,55].

In view of the fact there are no measurement of global direct or diffuse radiation at some stations in China, the 3-factor Re empirical model was adjusted to the 2-factor Re model, i.e., the scattering term was removed in Equation (2).

Secondly, nighttime Re, dark respiration was developed by further application of EMGPP (Equation (1)) or Re (LR) (Equation (2)) from daytime to nighttime and named the Re (DR) method.

E^−0.1bRetm × cos(Z) = C₁e^−kWm × cos(Z) + C₀

(3)

Similarly, we set b = 1, m = 1 and cos(Z) = 1 for nighttime Re simulation. Coefficients C₁ and C₀ were obtained by analyzing nighttime measurements (Re, water vapor). The data criteria as e^−0.1bRetm × cos(Z) > 1 was also used in nighttime Re model development. The hourly respiration during a day was computed by using above daytime respiration and nighttime respiration. This method is named Re (LDR). Similar methods are used in the studies [56,57,58].

Additionally, daytime Re was estimated by extrapolating nighttime EMRe (b = 1 and m = 1) and daytime parameters, which are widely used in current Re simulations [54,57]. We named this method as Re (DDR) for a whole day. In the following sections, evaluations of Re simulations in nighttime and daytime + nighttime using above methods were performed. So, an empirical model of respiration (EMRe) was used to simulate daytime and nighttime respiration.

Finally, NEP values in different time scales during 2013–2016 were calculated by GPP-Re and two Re models under all sky conditions.

3. Results

3.1. Model Development in a Subtropical Coniferous Forest

As PAR transfers in the atmosphere and its interactions with atmospheric GLPs are different under different atmospheric conditions, two situations were investigated for S/Q < 0.5 and S/Q ≥ 0.5, respectively, which represented a relative clean and a polluted atmosphere and environmental conditions (e.g., GLP concentrations, clouds, solar radiation, air temperature, relative humidity, and water vapor).

3.1.1. Empirical Model of Ecosystem Respiration in Daytime: Re (LR) Method

During the daytime and nighttime, Re and NEP displayed small variations in this coniferous forest compared to that of GPP. NEP showed evident diurnal, seasonal, and inter-annual variations, and was influenced by the driving factors, e.g., PAR, air and soil temperature, precipitation, clouds, water vapor, vapor pressure deficit (VPD), soil water content [46,54,59,60], and higher values appeared around noon and in the summer.

Selecting the hourly observational data described in Section 2.3 together with that corresponding to S/Q < 0.5 and S/Q ≥ 0.5, 67 and 985 group data (sample point n = 67 and 985) from 23 May 2013 to 31 December 2014 were used to develop the 3-factor and 2-factor Re empirical models during the daytime, respectively. Meanwhile, the water vapor pressure at the ground and Re were also used to determine the coefficients B_i of the EMRe models (defined as Re (LR)).

Table 1. Coefficients of the 3-factor and 2-factor EMRe models (3-F, 2-F) in the daytime (Equation (2)) determined using observational data in the subtropical coniferous plantation during 2013–2014, statistical metrics, i.e., coefficient of determination (R²), average of the absolute relative bias (δ_avg, %), normalized mean square error (NMSE), ratio of calculated Re to observed Re (cal/obs), standard deviations of calculated and observed fluxes (σ_cal and σ_obs), together with mean bias errors (MAD, mgCO₂ m⁻² s⁻¹ and %) and root mean square errors (RMSE, mgCO₂ m⁻² s⁻¹ and %) (S/Q < 0.5, n = 67).

Model	B1	B2	B3	B0	R2	δavg	NMSE	cal/obs	σcal	σobs	MAD		RMSE
											(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
3-F	0.026	1.036	0.000	0.057	0.999	24.31	0.085	1.001	0.07	0.02	0.05	24.04	0.06	29.37
2-F	0.026	1.036		0.057	0.999	24.31	0.085	1.001	0.07	0.02	0.05	24.04	0.06	29.37

and

Table 2. Same as Table 1, but for the 3-factor and 2-factor EMRe models at S/Q ≥ 0.5 (n = 985).

Model	B1	B2	B3	B0	R2	δavg	NMSE	cal/obs	σcal	σobs	MAD		RMSE
											(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
3-F	0.016	1.185	−0.112	0.093	0.994	68.45	0.582	1.009	0.17	0.04	0.12	63.41	0.14	76.72
2-F	0.008	1.234		0.039	0.994	72.46	0.638	1.010	0.18	0.04	0.12	66.35	0.15	80.30

show 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 Re values), normalized mean square error (NMSE = $\overline{{(y_{cal}-y_{obs})}^2}/({\overline{y}}_{cal}\times {\overline{y}}_{obs})$ [61]), 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 and Figure 2 show the scatter plots of the calculated versus the measured daytime hourly Re during 2013–2014 using the 3-factor and 2-factor Re models for S/Q < 0.5 and S/Q ≥ 0.5 situations.

In general, the 3-factor and 2-factor Re models displayed reasonable simulations for S/Q < 0.5 and S/Q ≥ 0.5 conditions, and larger estimation uncertainties were found for Re simulations at S/Q ≥ 0.5, which corresponds to large GLP loads, such as clouds, dust, aerosols, rain, dew, fog, resulting in large observational errors in flux, solar radiation, and meteorology measurements (Table 1 and Table 2). Under S/Q < 0.5 conditions (a much cleaner atmosphere compared to that at S/Q ≥ 0.5), the 3-factor Re model performed similar or the same estimations as the 2-factor model, indicating that the scattering factor/roles can be omitted.

3.1.2. Empirical Model of Ecosystem Respiration during Nighttime: Re (DR) Method

Using data criteria as e^−0.1bRetm × cos(Z) > 1 and in Section 2.3, 4412 and 4445 group data (n = 4412 and 4445) in 2013 and 2014 from 1 January to 31 December, along with the water vapor pressure at the ground and Re were selected to determine the coefficients C_i of the EMRe model (Equation (3)) during the nighttime. The results and statistical metrics are shown in

Table 3. Same as Table 1, but for the coefficients of the EMRe model (Equation (3)) determined using observational data (n = 4412, 4445) in 2013 and 2014, respectively.

Data	C1	C0	R2	δavg	NMSE	cal/obs	σcal	σobs	MAD		RMSE
									(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	0.153	−1.133	0.589	26.05	0.112	0.937	0.06	0.08	0.03	18.46	0.05	32.41
2014	0.175	−1.150	0.606	32.68	0.141	1.002	0.06	0.08	0.03	19.81	0.05	37.59

.

According to the better ratio of the calculated Re to the observed Re, 1.002, the coefficients of 2014 were selected for EMRe. Almost all the calculated Re values were close to the observed within one standard deviation, except that the observed Re values were smaller than −0.3 (absolute Re > 0.3). Figure 3 and Figure 4 display the calculated and observed hourly Re values in 2014 and their scatter plots. Generally, the EMRe simulations agreed reasonably with the observations during the nighttime in 2014.

It is noticed that there were some scattered points (Figure 4, left) that may influence the Re model’s performance. Therefore, seven larger scattered points in 2014 were removed in the further analysis, and the new estimated results are shown in

Table 4. Same as Table 3, but for using the updated data in 2014 (n = 4438).

Data	C1	C0	R2	δavg	NMSE	cal/obs	σcal	σobs	MAD		RMSE
									(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2014	0.175	−1.150	0.606	32.68	0.141	1.001	0.06	0.08	0.03	19.34	0.05	33.96

and Figure 4 (right). No large differences were found, e.g., coefficients and the ratio of the calculated Re to the observed Re. Then, the previous Re model was used in the later study.

3.2. Evaluation and Application of the EMRe Models under all sky Conditions during 2013–2016

3.2.1. Nighttime Re Simulations Using Re (DR) Method

We used the EMRe model (Equation (3)) for Re estimations from 22 May 2013 to 31 December 2016; the hourly Re values calculated using the coefficients in Table 3 were compared with the measurements. To fully evaluate the model performance, detailed results for each year and 2013–2016 are displayed in

Table 5. Same as Table 1, but for hourly nighttime Re during 2013–2016.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	27.24	0.121	0.064	0.074	0.92	2752	0.031	18.80	0.054	33.27
2014	44.76	0.148	0.065	0.083	1.01	4576	0.028	20.71	0.052	38.67
2015	44.40	0.110	0.057	0.070	0.93	4673	0.028	18.49	0.049	31.98
2016	43.46	0.227	0.060	0.095	0.81	4835	0.048	27.26	0.076	42.89
2013–2016	41.42	0.159	0.061	0.083	0.91	16,836	0.034	21.91	0.060	37.98

and Figure 5 (a scatted plot of calculated versus measured nighttime hourly Re during 2013–2016). Generally, hourly Re simulations exhibited suitable performances at nighttime with little underestimations for 2013 and 2015. A little larger underestimation (~20%) in 2016 may be caused by instrument problems during June–October.

Similarly, the daily sums of the nighttime Re values were also estimated and shown in

Table 6. Same as Table 5, but for daily sums of Re at nighttime during 2013–2016.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	15.35	0.029	0.656	0.553	0.92	224	0.248	12.34	0.329	16.37
2014	19.50	0.036	0.688	0.751	1.01	365	0.242	14.31	0.321	18.96
2015	12.14	0.024	0.602	0.582	0.93	365	0.210	10.73	0.290	14.78
2016	20.34	0.078	0.677	0.814	0.81	366	0.463	19.82	0.588	25.16
2013–2016	16.99	0.046	0.658	0.737	0.91	1320	0.296	14.79	0.407	20.36

. The daily sums of the Re simulations were also in line with the observations.

Furthermore, the mean monthly sums of the nighttime Re values for each year and four years over 2013–2016 were also computed (

Table 7. Same as Table 6, but for mean monthly sums of Re at nighttime during 2013–2016.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	11.69	0.013	45.22	41.61	0.92	7	5.517	9.02	7.041	11.51
2014	11.23	0.010	30.31	47.42	1.03	12	4.506	8.75	5.433	10.55
2015	7.75	0.008	41.05	41.06	0.93	12	4.244	7.12	5.221	8.76
2016	18.99	0.056	40.83	51.73	0.79	12	13.626	19.11	15.780	22.13
2013–2016	12.50	0.025	18.28	21.07	0.91	53	7.143	11.74	9.337	15.34

and Figure 6), and were also in agreement with the observations.

Figure 7 shows the simulations of the annual sums of Re at nighttime during 2013–2016. The ratios of the mean annual sums of the calculated to the observed Re were 0.93, 1.01, 0.93, and 0.81 for 2013, 2014, 2015, and 2016, as well as 0.91 for 2013–2016.

In a brief summary for the estimates of the nighttime Re from the hourly to annual time scales, the EMRe model exhibited a reasonable performance and most of the simulations were within a one time standard deviation of the observed Re.

3.2.2. Daytime and Nighttime Re Simulations Using Re (DDR) Method

Based on some studies where the daytime Re was extrapolated from the nighttime Re [37], the EMRe model at nighttime was further applied to calculate the Re at daytime during 22 May 2013–31 December 2016 using Equation (3) and observation data at daytime (water vapor). This Re (DDR) method was used for the Re simulation in the daytime and nighttime. The detailed results for the averages of the hourly, daily, monthly, and annual sums of Re are given in

Table 8. Same as Table 5, but for mean hourly daytime and nighttime Re during each year and 2013–2016.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	15.96	0.054	0.062	0.063	0.93	5376	0.021	12.46	0.038	22.35
2014	21.49	0.070	0.066	0.075	1.03	8755	0.022	15.32	0.039	26.90
2015	14.75	0.050	0.059	0.062	0.93	8753	0.020	12.24	0.035	31.52
2016	24.08	0.144	0.059	0.083	0.79	8722	0.044	23.25	0.064	33.73
2013–2016	19.40	0.159	0.062	0.074	0.91	31,606	0.027	16.46	0.046	27.71

,

Table 9. Same as Table 6, but for daily sums of Re at daytime and nighttime during 2013–2016.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	12.57	0.015	1.470	1.309	0.93	224	0.381	9.25	0.619	15.01
2014	15.71	0.018	1.570	1.555	1.03	365	0.369	10.74	0.608	17.72
2015	11.15	0.015	1.394	1.339	0.93	365	0.359	9.16	0.600	15.31
2016	21.99	0.077	1.398	1.681	0.79	366	0.945	21.03	0.974	21.66
2013–2016	15.66	0.035	1456	1.570	0.91	1320	0.528	13.27	0.708	17.70

and

Table 10. Same as Table 7, but for mean monthly sums of Re at daytime and nighttime during 2013–2016.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	10.64	0.009	45.22	41.61	0.92	7	9.952	7.95	7.041	11.51
2014	10.44	0.009	30.31	47.42	1.03	12	8.096	7.76	5.433	10.55
2015	8.06	0.007	41.05	41.06	0.93	12	8.419	7.06	5.221	8.76
2016	20.68	0.066	40.83	51.73	0.79	12	28.303	20.64	32.733	23.88
2013–2016	12.67	0.026	42.02	47.09	0.91	53	14.13	11.67	17.156	14.17

and Figure 8 and Figure 9.

For daytime + nighttime Re, (1) the EMRe also exhibited acceptable performances, and most of the estimations were within a one time standard deviation of the observed Re; (2) the estimation errors of Re were smaller than that that at nighttime (e.g., MAD, RMSE); thus, it may be caused by that σ_obs (observed Re) decreased with the increase in sampling number (i.e., from 16,836 to 31,606 for 2013–2016). Data quality of observed Re at daytime and nighttime, as well as only at nighttime in 2016, were the worst with the largest standard deviations (σ_obs, Table 5, Table 6, Table 7, Table 8, Table 9 and Table 10), which caused larger simulation errors in 2016 than in other years. The best Re estimates were performed for nighttime and daytime + nighttime Re in 2014 than in other years. It reveals that better observational data are very necessary for accurate Re evaluations.

For the consistency, the units used for the hourly, daily, monthly, and annual sums of Re were the same as mgCO₂ m⁻² s⁻¹. A coefficient can be multiplied by (3600 × 12/44)/1000 ≈ 0.982 to convert the hourly sums of Re to the widely used unit (gC m⁻²) for monthly and annual sums of Re.

Based on the above results, the Re values calculated by the EMRe model (Re (DDR) method) were multiplied by a factor 1.1 to best fit the hourly, daily, monthly, and annual Re values in averages and sums for each year and four years 2013–2016, so as to reduce the estimation error in the extrapolation of the respiration from nighttime to daytime; therefore, the updated Re estimates were improved at different time scales to some extent. As a brief summary, the ratios of the mean calculated to the observed Re were 1.000, 1.003, 1.001, and 1.001 for hourly, daily, monthly, and annual sums, respectively, from 2013 to2016. Then, the updated EMRe (i.e., using a modified coefficient 1.1) was applied in the later Re calculations.

3.2.3. Daytime and Nighttime Re Simulations Using Re (LDR) Method

As an alternative, daytime and nighttime Re values from 22 May 2013 to 31 December 2016 were computed using the 3-factor and 2-factor Re (LR) and Re (DR) methods (Equations (2) and (3)), respectively. This Re (LDR) method was applied to compute hourly Re in a whole day. The statistical metrics for the averages of the hourly sum of Re are shown in

Table 11. Same as Table 8, but for mean hourly daytime and nighttime Re during each year and 2013–2016 using 3-factor Re (LR). (Re (LDR) method).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013–2014	103.31	1.555	0.214	0.072	0.90	13,893	0.111	72.37	0.182	118.44
2015	77.37	1.165	0.202	0.062	0.97	8754	0.110	67.33	0.174	106.38
2016	78.60	1.401	0.208	0.083	0.73	8724	0.121	64.07	0.191	101.11
2013–2016	88.91	1.375	0.209	0.074	0.87	31,611	0.113	68.19	0.182	109.40

and

Table 12. Same as Table 11, but for using 2-factor Re (LR) (Re (LDR) method).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013–2014	104.24	1.530	0.210	0.072	0.86	13,893	0.106	69.30	0.176	114.70
2015	76.57	1.083	0.194	0.062	0.91	8754	0.103	63.09	0.162	99.32
2016	77.31	1.300	0.200	0.083	0.69	8724	0.112	59.17	0.179	94.93
2013–2016	88.65	1.309	0.202	0.074	0.83	31,611	0.107	64.17	0.173	103.93

and Figure 10. The ratios of the calculated to the observed Re values for the daily, monthly, and annual sums of Re were similar to that of the hourly sum (not given for saving space). Generally, the calculated Re values underestimated the observations by about 13% for the mean hourly and annual sums in 4 years.

3.2.4. Daytime and Nighttime Re Simulations Using Re (LDDR) Method

In view of the systemic underestimations of Re using the Re (LDR) method, dark respiration (Re (DR)) was considered in daytime Re simulation, i.e., using Re (LR) and Re (DR) methods to calculate daytime Re, respectively, and Re (DR) to calculate night Re. This method is named as Re (LDDR). The calculated results for the averages of the hourly sum of Re are displayed in

Table 13. Same as Table 11, but for mean hourly daytime and nighttime Re during each year and 2013–2016 using 3-factor Re (LDR) (Re (LDDR) method).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013–2014	112.13	1.586	0.256	0.072	1.40	13,891	0.144	93.82	0.223	149.16
2015	87.19	1.429	0.252	0.062	1.41	8753	0.141	86.01	0.232	141.95
2016	81.58	1.199	0.243	0.083	1.08	8722	0.137	72.86	0.214	113.53
2013–2016	96.61	1.414	0.251	0.074	1.31	31,606	0.141	85.06	0.226	135.86

and

Table 14. Same as Table 12, but for using 2-factor Re (LDR) (Re (LDDR) method).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013–2014	112.52	1.516	0.251	0.072	1.36	13,891	0.141	91.67	0.220	143.75
2015	86.10	1320	0.244	0.062	1.35	8753	0.136	83.33	0.218	133.48
2016	81.26	1.117	0.237	0.083	1.04	8722	0.134	70.87	0.204	107.93
2013–2016	96.46	1.333	0.245	0.074	1.26	31,606	0.138	82.92	0.216	129.77

and Figure 11. The ratios of the calculated to the observed Re values for daily, monthly, and annual sums of Re were also similar to that of hourly sum and not shown. In general, the simulated Re displayed evident seasonal variations, much higher in summer and lower in winter than the observations. The calculated Re values overestimated that measured by ~30% for mean hourly and annual sums in 4 years.

Several studies report the annual respiration biases of ~6% [62], 9.7% [31], and 20% [63,64], respectively. In general, the EMRe exhibited reasonable simulations of respiration.

3.3. NEP Simulations under All Sky Conditions during 2013–2016

The hourly net ecosystem productivity was calculated using the GPP-Re, EMGPP (3-factor and 2-factor, respectively, [45]) and EMRe models (i.e., Re (DDR), RE (LDR), Re (LDDR)), along with the observed PAR and water vapor pressure. To obtain better NEP estimates, several NEP calculation methods were used.

3.3.1. GPP Simulations Using EMGPP Model

Firstly, the hourly GPP estimates using the 3-factor and 2-factor EMGPP models are presented briefly in

Table 15. Same as Table 1, but for mean hourly GPP (mgCO₂ m⁻²) simulations (S/Q = 0–1) using 3-factor EMGPP during one year and 2013–2016.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	117.15	0.490	0.374	0.291	1.15	2493	0.287	59.43	0.535	111.05
2014	148.47	0.612	0.392	0.295	1.37	4181	0.301	71.90	0.549	130.99
2013–2016	142.74	0.586	0.381	0.300	1.31	14,748	0.318	70.03	0.398	87.60

and

Table 16. Same as Table 15, and for mean hourly GPP (mgCO₂ m⁻²) simulations (S/Q = 0–1) using 2-factor EMGPP.

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013	112.71	0.496	0.369	0.291	1.17	2493	0.281	58.31	0.530	110.00
2014	140.95	0.557	0.376	0.295	1.34	4181	0.289	69.00	0.536	127.94
2013–2016	135.55	0.538	0.366	0.300	1.29	14,748	0.303	66.62	0.378	83.23

[45]. Generally, the EMGPP model showed reasonable performance and overestimation of GPPs by ~30%.

3.3.2. NEP Simulations Using EMGPP and Re (DDR) Models

Secondly, the hourly NEPs were computed using EMGPP and Re (DDR) (with a modified factor 1.1). Generally, the ratios of the mean calculated to the measured Re in 2013, 2014, 2015, 2016, and 2013–2016 were 1.45, 1.89, 2.78, 3.12, and 2.25 for using the 3-factor EMNEP, and 1.54, 1.78, 2.60, 2.94, and 2.15 for using the 2-factor EMNEP. Similar ratios of the mean daily, monthly, and annual sums were also obtained. Thus, the NEP values were overestimated by a factor of 2.3 (ranged from 1.45 to 3.12).

3.3.3. NEP Simulations Using EMGPP and Re (LDR) Models

Thirdly, the hourly NEPs were calculated using EMGPP and Re (LDR) (with the 3-factor and 2-factor methods). The ratios of the mean calculated to the observed Re in 2013, 2014, 2015, 2016, and 2013–2016 were 1.94, 2.32, 2.97, 3.89, and 2.68 for using the 3-factor EMGPP and Re (LDR) models, and 2.04, 2.36, 3.01, 3.92, and 2.73 for using the 2-factor models. Thus, the hourly NEPs overestimated the observations by a factor of 2.7 (ranged from 1.94 to 3.92). Similar (but a little smaller) ratios of mean daily, monthly, and annual sums were also obtained.

3.3.4. NEP Simulations Using EMGPP and Re (LDDR) Models

Thirdly, the hourly NEPs were calculated using EMGPP and Re (LDDR) (with the 3-factor and 2-factor methods) and

Table 17. Same as Table 13, but for mean hourly NEP during each year and 2013–2016 (using 3-factor EMGPP and Re (LDDR) methods).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013–2014	755.05	9.008	0.245	0.293	0.93	13,894	0.108	174.46	0.178	288.95
2015	352.87	10.944	0.248	0.288	1.37	8760	0.106	238.23	0.173	387.05
2016	244.45	18.583	0.262	0.316	1.96	8777	0.131	392.98	0.202	604.30
2013–2016	497.35	11.281	0.251	0.299	1.22	31,671	0.114	230.34	0.183	371.58

and

Table 18. Same as Table 15, but for using 2-factor EMGPP and Re (LDDR) methods).

Time	δavg	NMSE	σcal	σobs	cal/obs	n	MAD		RMSE
Period							(mgCO2 m−2 s−1)	(%)	(mgCO2 m−2 s−1)	(%)
2013–2014	778.04	8.889	0.250	0.293	0.99	13,894	0.109	176.10	0.183	296.89
2015	361.45	10.906	0.251	0.288	1.41	8760	0.108	241.62	0.175	391.68
2016	246.72	19.059	0.263	0.316	1.97	8777	0.133	397.46	0.204	612.15
2013–2016	510.26	11.277	0.254	0.299	1.27	31,671	0.115	232.88	0.187	378.42

, and Figure 12 and Figure 13 gives the results. Similar ratios were found for the mean daily, monthly, and annual sums. The estimated NEPs in different time scales overestimated the observations by a factor of 1.22 (relative bias 22%) and 1.27 (relative bias 27%) for using the 3-factor and 2-factor models, respectively. The calculated standard deviation (0.251 mgCO₂ m⁻² s⁻¹) was less than that observed (0.299 mgCO₂ m⁻² s⁻¹) during 2013–2016 for using the 3-factor model, and the calculated standard deviation value was 0.254 mgCO₂ m⁻² s⁻¹ when using the 2-factor model.

It was found that the Re and GPP simulations captured the monthly variations well, but the NEP simulations did not show a better performance than Re and GPP, and further studies are needed. The estimates of the annual sums of NEP were in reasonable agreement with the observed NEPs, and most of them were within 2 σ_obs, except 2016.

Given there were problems in CO₂ flux measurements in 2016 summer, the observational data in 2013–2015 were used for further comparisons; the ratios of the mean annual sum of the calculated to observed values were 1.07 and 1.12 for NEP using the 3-factor and 2-factor methods, respectively, 1.41 and 1.36 for Re, and 1.07 and 1.12 for GPP. So, EMNEP models displayed improved performances for NEP. The NEP simulations exhibited the largest uncertainty than GPP and Re; this feature was also found in other studies, e.g., the Yale Interactive terrestrial Biosphere (YIBs model) [65], because NEP has the positive and negative values and the smallest values compared to GPP and Re, and the most important reason is that NEP is associated with many complex processes and current limited understanding.

According to the high NEP in the summer of 2013–2015, the measured low NEP in 2016 summer should be corrected, and the simulated NEP values can be used as a reference.

4. Discussion

4.1. PAR Balance Principle

The emissions of the biogenic volatile organic compounds (BVOCs) from vegetation are part of carbon emissions, and play important roles in the atmosphere, biosphere, and climate change [66,67,68,69,70,71,72,73,74,75,76,77]. Similar to the CO₂ flux, BVOC emissions are controlled by key driving factors, PAR, air/leaf temperature, water, etc., and are suggested to be considered in carbon balance.

EMRe expresses PAR balance and usage above a canopy level, associating with the PAR attenuation of the equivalent CO₂ for respiration, GLPs’ (except CO₂) absorption and use of PAR, and GLPs scattering of PAR. Based on the above Re and NEP simulations in different time scales, including Re simulations in daytime and nighttime, further applications of PAR balance and EMBE were acceptable. In summary, PAR energy and its use play significant roles in the processes of CO₂ and BVOCs, i.e., PAR controls the processes of GPP, Re, and NEP in different ways/contributions, which are exhibited and represented by each coefficient and constant in the empirical models/energy relationships. Therefore, PAR provides a basic and vital energy for all biogenic activities and processes. Furthermore, PAR simultaneously interacts with all physical, chemical, and biological processes in the atmosphere-vegetation-land system, i.e., not in only one direction but in two directions and several directions among the processes in the atmosphere, vegetation, and land [45,55]. Namely, direct and/or indirect PAR use is applicable to the study of BVOCs, GPP, Re, and NEP, and one basic principle of PAR utilization runs through them all.

PAR energy method is a universal and practical method and deserved to be investigated in the future. It describes main processes through energy and saves time to deal with numerous specific processes in the atmosphere-vegetation-land system. The advantage is that this method can reveal multiple interactions and the mechanisms associated with CO₂ and BVOCs and their main driving factors (temperature, PAR, water, scattering of solar radiation S/Q). As an example, the universality of PAR energy use/method was also exhibited in the Re empirical model development that extrapolated form daytime to nighttime (Equations (2) and (3), Section 3).

4.2. The Application of Re (DR) Model and Calibration of Longwave Radiation Sensor

The PAR balance equation can be obtained from EMGPP (Equation (1)) by replacing GPP with isoprene (one dominant species of BVOCs), or converted from the EMBE (Section 5.4 in [55]):

PAR = C₁e^−a1^ISOm × cos(Z) + C₂e^−kWm × cos(Z) + C₃e^−S/Q + C₀

(4)

The computed PAR at the top of the atmosphere (TOA) is C₁+ C₂+ C₃ + |C₀| and is equal to the PAR at the TOA measured by the satellites and airplanes (PAR_TOA). This method can be used to calibrate the PAR sensor by the ratio of PAR_TOA/(C₁+ C₂+ C₃ + |C₀|).

Similarly, the Re (DR) model (Equation (3)) also exhibited an energy balance at the canopy level. Based on the Beer–Lambert law, C₁ represents the longwave radiation at the TOA associated with the individual atmospheric GLPs, and C₀ represents that emitted from the surface (i.e., vegetation). The processes and their associated energy are divided and quantified by separate terms in empirical equations. So, C₁ +|C₀| = 357.22 and 368.06 W m⁻² is the estimated total longwave radiation contributed by the atmosphere and the surface for 2013 and 2014 night (using half-hour flux and other corresponding data in this section, n = 4412, 4445, respectively). The detailed meaning of C₁ +|C₀| is described in studies [55,78]. It should be noted that the longwave radiation of the GLPs at the TOA is isotropy. According to the canopy temperature measured at the flux tower (23.6 m), forest emissivity = 0.95 [79] and Stefan–Boltzmann law, the estimated longwave radiation emitted by vegetation canopy were 379.87 and 379.91 W m⁻² with relative biases of 5.96% and 3.12% for 2013 and 2014 night, respectively. This estimated longwave radiation can be used as a reference at TOA, and the ratios of 0.940 and 0.969 can be selected as calibration factors. Compared this estimated total longwave radiation to that emitted from the canopy measured at the tower (41.6 m) (398.19 and 399.29 W m⁻²), their relative biases were 10.28% and 7.82%, respectively.

The attenuating rates of CO₂ and photochemical terms using Equation (3) were 1.39% and 77.60% for 2013, respectively, and 1.42% and 77.40% for 2014, indicating the longwave radiation absorption by CO₂ was much lower than that of atmospheric components, and can be neglected, considering its low concentration compared with other GLPs. These were in good agreement with the studies [80,81].

The annual averages of downward longwave radiation at the canopy level were computed using Equation (3), and a longwave radiation that was remitted by the atmospheric GLPs after they absorb the longwave radiation emitted from the canopy, |C₀| × (1 − e^−kWm × cos(Z)) were added. Then, their annual averages were 352.30 and 358.55 Wm⁻² for the 2013 and 2014 night, respectively, the corresponding values measured at 41.6 m were 362.31 and 362.01 Wm⁻², respectively. Their relative biases were 0.29% and 1.22% for 2013 and 2014, respectively.

In a summary, the Re (DR) model showed a good performance for the estimations of the downward and upward longwave radiation. It was a new piece of evidence to demonstrate the reasonability and potential application of the Re (DR) model.

4.3. EMNEP Model and Its Simulation

The hourly NEPs were calculated by GPP-Re using the 3-factor and 2-factor EMGPP and EMRe. It is clear that the NEP performance in the hourly, monthly, and annual NEP simulations were gradually improved when more processes in respiration were considered (e.g., the bias and RMSE values were decreased, Section 3.3). The purpose of developing four types of Re empirical models (Section 3.2) was for the better simulation of NEP, one of the most important factors in carbon exchange, and the best NEP estimations, were found by using Re (LDDR) and the corresponding GPP models (Section 3.3). Finally, the NEP empirical model captured most of the hourly, monthly, and annual NEP variational patterns (within 2 σ_obs). The CO₂ flux measurements in 2016 (July, August, and October) had some problems and caused many unusual drops in the monthly sums that exceeded 2 σ_obs; therefore, the simulated NEPs can be used as a reliable reference for these problematical periods in the CO₂ flux measurements. Generally, the monthly sums of NEP exhibited overestimates in the summer and underestimates in the winter between 2013 and 2016, and the annual sums of NEP were in 1 σ_obs for 2013 and 2014 and 2 σ_obs for 2015.

Under S/Q <0.5 conditions (low GLP loads, low cloud amounts, and aerosol concentrations), the constant in the Re (LR) method was close to 0, indicating that the PAR attenuated by atmospheric substances plays absolute dominant roles in respiration, while the scattered PAR does not. Under S/Q ≥ 0.5 conditions (high GLP loads, cloud amounts, and aerosol concentrations), both the attenuated and scattered PAR by GLPs contribute respiration, or both the direct and diffuse PAR contribute the respiration and the diffuse PAR contributes much more than that at S/Q <0.5. A similar result is also reported [82,83,84].

The respiration is small compared to NEE and GPP, and very difficult to be simulated accurately. The improved performance in the NEP simulations (hourly to annual sums) from the Re (LR) to Re (LDDR) method revealed that some factors (dark respiration and photorespiration, etc.) should be taken into consideration in the respiration and NEP processes. Based on these Re calculation results, the different NEP estimates can be used as references to evaluate NEP uncertainties to some extent.

4.4. GPP, Re, and NEP Empirical Models

Re and NEP empirical models captured some basic characteristics of Re and NEP, which were in line with other studies. It should be noted that GPP, Re, and NEP empirical models describe a dynamic and multiple-direction energy balance, more specifically, quantify the relationships between energy (PAR) and atmospheric substances (CO₂ and other absorbing and scattering GLPs). The complicated and interactive CO₂ processes in the atmosphere, biosphere, and land were expressed using PAR and its related energy. It is an advantage that these PAR energy empirical models can simplify numerous CO₂ processes and capture the main variation properties of CO₂ (GPP, Re, and NEP). Thus, we do not need to pay too much attention to specific processes. In the future, to obtain more reasonable simulations of seasonal NEP, more studies are necessary. It should be emphasized that solar radiation, atmosphere, biosphere, and land should be considered as a whole system. The interactions in PAR energy usage between photons, CO₂ molecules, and other GLPs quantum obey the statistical law. It should be noted that there were larger RMSE values for the simulations of Re in 2016 (e.g., Table 8, Table 9 and Table 10) than in other time periods, which was also mainly caused by measurement problems and very low Re. Thus, timely instrument maintenance is significantly important, and the empirical model can provide valuable data.

Water vapor is an objective parameter to represent air temperature and humidity (T, RH), as water vapor is a function of T and RH [85]. Water vapor is also a vital participant and reactor in chemical and photochemical reactions, associating with OH radicals and BVOCs, as well as solar radiation [55,86].

The simulated monthly sums of Re and NEP exhibited larger variation ranges than their observations (Section 3.2), which may reveal some problems in data processing methods, e.g., spike detection, threshold (u*_c) determination of nocturnal friction velocity (u*), and gap-filling model selection had important impacts on the estimation of annual NEE [87]. In more detail, the relative estimation deviations of annual gross ecosystem exchange and ecosystem respiration induced by the uncertainty in data processing were 3.88–11.41% and 6.45–24.91%, respectively [87]. Additionally, more than 50% of the night data are eliminated using the threshold (u*_c) [87]. Large original flux data (e.g., the assumed abnormal data) removal caused by different methods may result in small variation ranges in measured Re and NEP and inaccurate understanding. The ratios of the observed annual sum Re/GPP and NEE/GPP were 0.764 and 0.236, respectively, for 2013–2016. These small ratios help us to clearly understand the different levels and ranges of GPP, Re, and NEP. Therefore, more and more strict appreciated data processing methods are required to improve flux data quality.

4.5. Some Issues Related to the Empirical Models of GPP, Re and NEP

Briefly, PAR absorption and indirect use by atmospheric components through chemical and photochemical reactions mainly associated with excited NO₂ (NO₂*), H₂O, OH radicals [88] and BVOCs + OH + O₃. Thus, the photochemical term expressed the total PAR used by GLPs. In Equations (1) and (2), the multiple interactions between PAR and the terms (photochemical, scattering, Re term in this study instead of BVOC terms) are explained in EMBE development in more detail [55,89].

The EMGPP showed that the GPP increases with the decrease in S/Q (i.e., the increase in the scattered PAR by GLPs, including aerosols and clouds) [45], which were also in good agreement with model studies and observations (the enhanced diffuse radiation caused by the aerosols results in the enhancement of GPP) [84,90,91,92,93,94,95,96,97]. Considering similar phenomena for respiration (Section 3.1.1) together, the diffuse radiation associated with the atmospheric GLP substances, i.e., light use efficiency, how much direct and diffuse PAR is absorbed by plants, play significant roles in plant growth and model simulations.

To better evaluate the empirical model performance, it is beneficial to understand the uncertainties in the measurements and popular model predications. A relative uncertainty of 20–50% for the CO₂ eddy covariance measurements is reported in [29]. The annual sums of NEE are 600 ± 130 gC m⁻² year⁻¹, 300 ± 180 gC m⁻² year⁻¹ at different forests [26]. Uncertainty for the annual NEE estimation is 12–32% [27]. A range of ±32% for the annual NEE, GPP, and ecosystem respiration values was used [98]. The annual NEP estimates using traditional ecological models agree with the flux measurements in 5%, 30%, and 100% [24,25,26]. The relative error of the PAR sensor is better than ±5% [78]. The yearly non-stability is <2% for global solar radiation sensor and <1% for solar direct radiation sensor [86]. The averaged relative bias of water vapor is less than 18% [99]. The uncertainties/responses of GPP to its driving factors (PAR, water vapor, and S/Q) using EMGPP is fully reported in [45] and can be used as a reference for the Re and NEP simulations.

Using the PAR energy balance method in this coniferous forest, we made applications of EMBE to GPP, then the Re and NEP simulations. Generally, the reasonable estimates of Re and NEP were obtained. The GPP, Re, and NEP empirical models are simple, timesaving, and the input parameters are easily obtained from the weather and solar radiation stations. They also revealed some basic mechanisms, which agreed with model studies, e.g., GPP, Re, NEP, and NEE increased with the increases in the diffuse PAR [82,83]. It is an important application that the Re (DR) model can also be used to calibrate the longwave radiation sensor (Section 4.2). In all, the empirical models in simulations of BVOCs, GPP, Re, and NEP are useful and suggested to be investigated at other flux stations for a more thorough estimation of the total carbon balance [70,100].

5. Conclusions

Application of PAR balance at the canopy level, 3-factor and 2-factor Re empirical models to describe respirations at different situations were developed at a subtropical coniferous plantation in China. Reasonable Re simulations were obtained in different time scales (hourly, daily, monthly, and annual) during 2013–2016. Then, hourly, daily, monthly, and annual sums of NEPs were calculated using the combined empirical models of GPP and Re, and the simulated NEPs also reasonably agreed with the measurements (within 2 σ_obs). The Re and NEP empirical models captured most variation patterns of Re and NEP in 2013–2016 and have better performances for annual sums of Re and NEP In general. As an example, in 4 years, the simulated Re values overestimated the observations by about 30% for mean hourly and annual sums. Similarly, the computed NEPs also overestimated the observations by about 22% and 27% in different time scales for using the 3-factor and 2-factor models, respectively, combined using previously developed 3-factor and 2-factor GPP empirical models. Empirical model simulations of RE, NEP, and GPP can be used to fill the gaps of RE, NEP, and GPP caused by the problems in CO₂ flux measurements. The Re (DR) model has another application for the calibration of longwave radiation sensors.

The PAR energy balance well-displayed a “universal” principle for expressing the energy relationships/main processes associated with Re and NEP, as well as GPP and BVOCs. It can simplify the complex processes of the CO₂ in the atmosphere and biosphere and estimate our concerned key CO₂ parameters. The PAR energy balance method exhibited some advantages, revealed some phenomena and mechanisms associating with Re and GPP that were in line with other extensively used models’ studies. It helps us to understand the complicated and multiple interactions/processes fully and deeply among solar radiation, atmosphere, biosphere, and land. Extensive studies deserved to be carried out at this site and other sites in the future.