Assessment of SITE for CO₂ and Energy Fluxes Simulations in a Seasonally Dry Tropical Forest (Caatinga Ecosystem)

Abstract

Although seasonally dry tropical forests are considered invaluable to a greater understanding of global carbon fluxes, they remain as one of the ecosystems with the fewest observations. In this context, ecological and ecosystem models can be used as alternative methods to answer questions related to the interactions between the biosphere and the atmosphere in dry forests. The objective of this study was to calibrate the simple tropical ecosystem model (SITE) and evaluate its performance in characterizing the annual and seasonal behavior of the energy and carbon fluxes in a preserved fragment of the Caatinga biome. The SITE model exhibited reasonable applicability to simulate variations in CO₂ and energy fluxes (r > 0.7). Results showed that the calibrated set of vegetation parameters adequately simulated gross primary productivity (GPP) and net ecosystem CO₂ exchange (NEE). The SITE model was also able to accurately retrieve the time at which daily GPP and NEE peaked. The model was able to simulate the partition of the available energy into sensible and latent heat fluxes and soil heat flux when the calibrated parameters were used. Therefore, changes in the dynamics of dry forests should be taken into consideration in the modeling of ecosystem carbon balances.

1. Introduction

With the intensification of global climate change, phenomena such as El Niño and La Niña are becoming more frequent and progressively affecting larger areas [1]. Simulations show that global climate change can adversely affect arid and semiarid regions by extending the dry period and jeopardizing the biodiversity of the ecosystem [2,3].

Seasonally dry tropical forests (SDTF) occur widely in the world, mainly in the Americas, Central and South Africa and India. The semiarid lands of Brazil are widely dispersed throughout the North-Eastern region, covering an area of approximately 865,000 km² with the Caatinga biome as the main vegetation type [4,5]. The Caatinga is a biodiversity-rich area [6], but its landscape is changing due to intense human activities, particularly deforestation and fires (slash-and-burn practices), with large native vegetation areas being replaced by pastures [5]. Such a process has important implications in land-use changes, and consequently in mass and heat fluxes changes in the soil–vegetation–atmosphere interface [7]. Thus, it is relevant to understand and quantify processes related to the Caatinga energy balance and carbon flux, which can be essential for the formulation of environmental and public policies.

Forests are widely considered to have the largest potential to act as sinks for atmospheric CO₂ [8]. Several studies have shown the potential of different types of ecosystems to sequester carbon in response to the increased atmospheric CO₂ associated with climate change or to release CO₂ as a result of changes in land use and management [9,10,11,12]. In Brazil, energy and CO₂ fluxes have been studied in biomes such as wetlands—the Pantanal [13]; Cerrado [14]; the Caatinga [15,16], but mostly in the Amazon region [17,18,19,20].

The environmental biology of the Caatinga vegetation (phenology, energetic seasonality, carbon and water fluxes) is a key component in the study of regional carbon cycle [4,16,21,22] and for the understanding of the ecosystem dynamics and how environmental drivers may affect it (e.g., recovery after wildfires and land use changes).The importance of the Caatinga for the regional climate [23,24] is therefore unquestionable [5,25], but its vulnerability to drought and the risks associated with a drier climate [26] is uncertain, due to a lack of studies in this region.

The eddy covariance technique (EC) has been widely used to determine mass and heat fluxes in the soil–vegetation–atmosphere interface. It is particularly useful to assess whether a particular ecosystem is acting as sources or sinks for atmospheric CO₂ [27,28,29]. Flux measurements also provide data for the evaluation of dynamic climate models [30]. EC is also useful to compare soil–vegetation–atmosphere transfer (SVATS models) with the observed seasonal behavior of ecosystem–atmosphere exchanges [31].

Few studies on the partitioning of the energy balance and/or carbon dioxide and water fluxes have been conducted in the Caatinga, particularly in areas of preserved or in-recovery vegetation [15,16,32,33,34], due to the complexity and cost of the field experiment installations. In this context, ecological and ecosystem models can be used as alternative methods to answer questions related to the interactions between the biosphere and the atmosphere. However, it is crucial to compare and validate the simulated data with those measured in situ in order to evaluate the simulations accuracy, and to calibrate specific parameters of the model according to the characteristics of a given ecosystem. The main challenge for the modeling of CO₂ and energy flux in the Caatinga is the adjustments of the site-specific biophysical and morphophysiological vegetation characteristics inherent to the different seasons of the year. For example, increasing the leaf area index (LAI), specific leaf area (SLA) and dimension of leaves may increase the rate of photosynthetic CO₂ fixation [9]. However, that increase in LAI and SLA can also overestimate water loss through evapotranspiration and, consequentially the latent heat flux may increase especially in the seasons with low water availability.

The objective of this study was to calibrate the simple tropical ecosystem model (SITE) and evaluate its performance in characterizing the annual and seasonal behavior of the energy and carbon fluxes in a preserved fragment of the Caatinga biome. The simple tropical ecosystem model (SITE) was used because it is a biogeochemical and biophysical model developed to simulate mass and energy fluxes between the ecosystem and the atmosphere, adopting ecosystem equations integrated over time. The SITE model is considered to be of intermediate complexity, sophisticated enough to be used in the study of the fast dynamics of tropical ecosystems [35]. The advantage of this model is that it can realistically represent complex interactions between precipitation, horizontal wind speed, light, temperature, and air humidity by incorporating observational data as input variables.

2. Materials and Methods

2.1. Description of the Experimental Area

The study was conducted in a preserved fragment of a seasonally dry tropical forest, the Caatinga biome, in the semiarid lands of the Northeast Brazil (6°34′42″ S, 37°15′05″ W, 205 m above sea level). The vegetation is composed by deciduous and semi-deciduous species with a shrub-tree structure, with approximately 8 m of height, predominantly sparsely distributed small trees and shrubs besides herb patches which thrive only during the wet season [36]. The region has a mean air temperature of 25 °C with little variability throughout the year. Mean annual rainfall ranges from 300 mm to 1000 mm, concentrated mostly in a 3–5 months period (January to May), followed by an extended dry season lasting from 7 to 8 months (June to December), with an average annual relative humidity of the air around 60 % [37,38]. Evapotranspiration rates are high (between 1500 and 2000 mm year⁻¹) and the predominant soil type is a sandy loam and sandy clay loam Entisol, shallow, rocky, with low fertility and low water holding capacity [38,39].

2.2. Micrometeorological Measurements

The data used for the calibration and validation of SITE were obtained through a micrometeorological tower of 11 m height, equipped with an EC system installed in a conservation unit of the Caatinga biome named Seridó Ecological Station (ESEC-Seridó), near the town of Serra Negra do Norte, in the Rio Grande do Norte State. The ESEC-Seridó area is managed by the Chico Mendes Institute for Biodiversity Conservation (ICMBio), comprising an area of 1163 ha of remnant Caatinga, which can be considered representative of the whole Caatinga biome. The micrometeorological tower belongs to the Brazilian National Institute of Semiarid (INSA) and is part of the National Observatory of Water and Carbon Dynamics in the Caatinga Biome (NOWCDCB) network. The data were collected from January to December 2014 and 2015.

In the study region, annual accumulated rainfall in 2014 and 2015 was of 513 mm and 466 mm, respectively while the 30-year average in the region is of 758 mm. Highest rainfall amounts were observed during the months from March to May and the highest daily value was 42 mm in 2014 and 65 mm in 2015, which clearly indicates the occurrence of a more intense and longer dry season. Mean annual temperature was 28.9 °C, ranging from 24.7 °C in March to 32.2 °C in November. Mean soil temperature was 31.4 °C, and annual integrated Rg was 8030 MJ m⁻². The mean annual value for the VPD was 1.7 kPa. The highest daily VPD means were registered in the dry season (2.7 kPa) while the minimum means (0.2 kPa) were registered in the wet season.

The measured data were grouped into two data sets: high and low frequency. The high frequency data obtained through the EC method consist of measurements of CO₂ and water vapor concentration and the three wind components ($u_x$, $u_y$, $u_z$) measured by an Integrated CO₂/H₂O Open-Path Gas Analyzer and 3D Sonic Anemometer (IRGASON, Campbell Scientific, Inc., Logan, UT, USA). Atmospheric pressure was measured using an Enhanced Barometer PTB110 (Vaisala Corporation, Helsinki, Finland) and air temperature was measured by a HMP155A probe (Vaisala Corporation, Helsinki, Finland). All these measurements were collected and stored at 10 Hz frequency in a memory card attached to a CR1000 model Datalogger (Campbell Scientific, Inc., Logan, UT, USA).

The CNR4 Net-radiometer (Kipp and Zonen, Delft, The Netherlands) was used to obtain low frequency data such as net radiation (incoming and outgoing shortwave radiation, longwave radiation emitted and reflected by the surface and longwave radiation emitted by the atmosphere) and albedo. Air temperature and relative humidity were measured with a HMP45C probe (Vaisala Corporation, Helsinki, Finland). All sensors were installed at a height of 11 m above the surface, around 4.0 m above the average vegetation canopy of the region. Soil heat flux density (G) was obtained through the average value between the measurements of two HFP01SC model plates (Hukseflux Thermal Sensors, Delft, The Netherlands), installed at a depth of 0.05 m. All data were sampled at a 5 s frequency and stored as half-hourly averages.

2.3. Data Processing and Post Processing

The fluxes of energy and CO₂ were calculated using the LoggerNet software (Campbell Scientific, Inc., Logan, UT, USA) by converting the high frequency data into the binary format (TOB1) with a 30 min timestep. The high frequency data were processed using the EdiRe software, developed by John Moncrieff and Robert Clement of the University of Edinburgh. The EdiRe algorithm transforms high frequency data in half-hourly averages. Data processing includes corrections such as: the detection of spikes, delay correction of H₂O/CO₂ in relation to the wind vertical component, coordinates rotation correction (2D rotation) using the planar fit method, correction of spectral loss, sonic virtual temperature correction, corrections for flux density fluctuation: WPL-correction [40], as well as the incorporated frequency response correction derived from the following studies [41,42].

For the detection of spurious data (spikes) we used a method described in the literature [43]. Gap filling due to data inconsistency and the rejection of spurious values was carried out by using the method described as proposed in the literature [44], which takes into consideration the covariance between fluxes and meteorological variables and also the temporal self-correlation of fluxes. In this algorithm, actions are taken considering the following conditions: (i) if there are missing flux data, but meteorological data (incoming solar radiation—Rg, Ta and vapor pressure deficit—VPD) are available, then the gap is filled with the mean value considering similar meteorological conditions in a 7-day window; (ii) if only incoming solar radiation data are available, the gap is filled with the mean value considering similar meteorological conditions in a 7-day window; (iii) if no meteorological data are available, the gap is filled with the mean value in the last hour, thus considering the diurnal variation of each variable. If data gaps still exist after applying the algorithm, the same procedures will be carried out but considering larger time windows. The gap filling method was carried out by using an online tool by the Max Planck Institute. Further details about data processing and post processing were presented in the literature [15].

2.4. Energy Balance

The energy balance equation expresses the conversion of net radiation (Rn) into energy and mass fluxes between the surface and the atmosphere:

$$\begin{array}{ccc}Rn=LE+H+G& \to & \left(W m^{-2}\right)\end{array}$$

(1)

where $LE$ is the latent heat flux density, $H$ is the sensible heat flux density and $G$ is the soil heat flux density, respectively, expressed in W m⁻². Sensible and latent heat turbulent fluxes were determined using high frequency data measured by the eddy covariance system based on the equations:

$$\begin{array}{ccc}LE={\rho }_{air}ƃ\lambda ·\overline{w^{\prime }{q}^{\prime }}& \to & \left(W m^{-2}\right)\end{array}$$

(2)

$$\begin{array}{ccc}H={\rho }_{air}·c_P·\overline{w^{\prime }{T}^{\prime }}& \to & \left(W m^{-2}\right)\end{array}$$

(3)

where ${\rho }_{air}$ is the density of the air, $\lambda$ is the latent heat of water vaporization, $c_P$ is the specific heat at constant pressure, $\overline{w^{\prime }{q}^{\prime }}$ and $\overline{w^{\prime }{T}^{\prime }}$ are the covariances between the deviations in vertical wind speed $\left(w^{\prime }\right)$ and the deviations in specific humidity $\left(q^{\prime }\right)$ and air temperature $\left(q^{\prime }\right)$.

2.5. Net Ecosystem Exchange

Net ecosystem exchange (hereafter referred to as NEE) is the sum of the eddy CO₂ flux $\left(F_{C{O}_2}\right)$ calculated as the covariance between the fluctuations of the vertical wind speed $\left(w^{\prime }\right)$ and the density of CO₂ $\left(c^{\prime }\right)$, and the rate of change of CO₂ stored in the air column below the EC measurement height $\left(Sc\right)$. Since no concentration profile was installed at the site, we opted for the discrete approach, assuming that the CO₂ concentration inside the canopy can be estimated as an approximation [45]. Thus, the half-hour values of the $Sc$ were calculated using the method proposed by [46] and widely used in subsequent studies [45,47].

The CO₂ fluxes were partitioned in order to separate NEE into GPP and ecosystem respiration (R_eco). We used a method of flux partitioning based on night-time values as proposed by [44]. For nocturnal periods, we considered the GPP equal to zero and therefore the NEE was estimated as follows:

$$NEE=R_{eco} \to \mathrm{for} \mathrm{night}-\mathrm{time} \mathrm{periods}$$

(4)

$$NEE=R_{eco}-GPP \to \mathrm{for} \mathrm{daytime} \mathrm{periods}$$

(5)

Night-time fluxes were adjusted according to air temperature (T_a) using the Lloyd and Taylor equation [48]:

$$R_{eco}=R_{eco.ref}*e^{E_0\cdot \left(\frac{1}{T_{ref}-T_0}-\frac{1}{T_a-T_0}\right)}$$

(6)

where, R_eco (μmol m ⁻² s ⁻²) is the sum of autotrophic and heterotrophic respiration, R_eco._ref is the respiration rate at a reference temperature T_ref (15 °C), E₀ (K) is the activation energy or the dependence of R_eco on temperature expressed in a temperature scale and T₀ is the base temperature set to −42.02 °C as suggested in the literature [48]. This model relates R_eco to T_a for night-time data and the temperature response function is then used to extrapolate R_eco values for the daytime periods. R_eco and GPP were calculated using the online tool provided by the Max Plank Institute for Biogeochemistry [49].

2.6. Description of SITE Model and Site Specific Biophysical Parameters

The 30 min micrometeorological observations measured at the study site were used to drive SITE simulations (SITE, version 1.1–0d). Observations included air temperature, precipitation, horizontal wind speed, downward shortwave and longwave radiation and specific humidity. The model was developed by Santos and Costa (2004), and is available at [50].

SITE is available as a FORTRAN code and structured with a canopy layer and two layers of soil using data collected above the canopy in 1 h periods, such as air temperature and specific humidity, horizontal wind speed, incident shortwave radiation, net longwave radiation, albedo, precipitation and atmospheric pressure. The model also uses parameters of biophysical characteristics of the vegetation. The main output variables of the model are net radiation (Rn, W m⁻²), latent heat flux (H, W m⁻²), sensible heat flux (LE, W m⁻²), soil heat flux (G, W m⁻²), gross primary production (GPP, g C m⁻² h⁻¹) and net ecosystem exchanges (NEE, kg C ha⁻¹ h⁻¹).

SITE is a dynamical point model that uses an integration time step (dt) of one hour, representing a land portion entirely covered by an evergreen broadleaf forest. The biophysical parameters used for the studied area were determined by in situ analyses and previous studies carried out in the Caatinga vegetation, shown in

Table 1. Site-specific parameters used in the simulation.

Parameter	Used Value	Source
Height of data measurement (z)	11 m	Measured on site
Height of the canopy (z1)	8 m	Measured on site
Height of lower canopy (z2)	5 m	Measured on site
Zero plane displacement (d)	7.33 m	Estimated
Roughness above the canopy (zh)	1.35 m	Estimated
Total soil porosity (Φ)	0.41 m3 m−3	[38,39]
Humidity content at field capacity (θCC)	0.225 m3 m−3	[38,39]
Moisture content of the permanent wilting point (θPM)	0.151 m3 m−3	[38,39]

. The SITE model uses a parameterization for the carbon balance adapted from the integrated biosphere simulator—IBIS model [7]. The hydraulic parameters of the soil were calculated as adapted from [51]. For more details on the dynamics of the model see the study [35].

A series of 730 days of data, collected from January to December 2014 and 2015, was used to evaluate the model. Moreover, a series of 120 days of data collected during the wet season (February to May 2014) and 92 days during the dry season (August to October 2014) was used to calibrate daily variations in the model. Validation of the model was performed using a dataset measured by the EC technique. The measured data used for calibration and validation of the model were: net radiation (Rn), latent heat (H), sensible heat (LE), soil heat (G) and CO₂ fluxes.

The SITE model was used using site-specific data on temperature, horizontal wind speed, specific humidity, photosynthetically active radiation (PAR) and rainfall in 2014 and 2015. The daily cycle and seasonal dynamics of simulated Rn, LE, H, G, GPP and NEE from the SITE model were compared with observed data form the wet and dry seasons. The SITE model was originally developed to study the response of tropical ecosystems to varying environmental conditions. Since our study area is a seasonally dry tropical forest with a long dry season, we carried out calibrations tests to assess the variability in energy flux and CO₂ caused by the variation in biophysical parameters.

Based on previous sensitivity analysis [52,53,54], we chose to calibrate only the parameters most likely to influence the result of the model: initial fraction of moisture in the soil (θg/θd), coefficient of stomatal conductance (m), maximum capacity of the Rubisco enzyme (V_max), typical dimension of leaves (d_u), typical dimension of stems (d_s), leaf width (w), and specific leaf area (SLA). Parametrization of the model followed the original parametrization included in the model SITE. Table 1 shows the parameters adopted for the Caatinga biome.

For the calibration of the model parameters we adopted the sequential method, in which parameters are calibrated separately for each step, according to the hierarchy used by the SITE model in the calculations: infrared radiation balance in the canopy and balance of solar radiation, aerodynamic processes, plant physiology, transpiration, balance of water intercepted by the canopy, transport of mass and energy fluxes, soil heat flux and soil moisture and carbon balance. In the sequential calibration procedure, we varied parameters in each step individually while all other parameters remained unchanged, as previously carried out by [55,56].

After exhaustive calibration tests of the SITE model for the Caatinga vegetation, we defined the confidence interval of the calibrated parameters, which is crucial to obtain reliable energy and carbon fluxes values. Afterwards, we ran 44 simulations varying the most relevant calibration parameters for the wet season and dry season, as follows: initial fraction of moisture in the soil (θg/θd) from 0.05 to 0.36, coefficient of stomatal conductance (m) from 4 to 10, maximum capacity of the Rubisco enzyme (V_max) from 40 to 120 μmol CO₂ m⁻² s⁻¹, typical dimension of leaves (d_u) from 0.020 to 0.090, typical dimension of stems (d_s) from 0.010 to 0.15, leaf width (w) from 0.01 to 0.15, and specific leaf area (SLA) from 9 to 26 m² leaf kg⁻¹ C. The ranges of parameter values were chosen according to either in situ measurements or previously published results in the Caatinga biome as described in

Table 2. Calibrated parameters: Specific leaf area (sla, m² leaf kg⁻¹ C), typical dimension of leaves (d_u, m), typical dimension of stems (d_s, m), leaf width (w, m), coefficient of stomatal conductance (m, dimensionless), maximum capacity of the Rubisco enzyme (V_max, μmol CO₂ m⁻² s⁻¹), and initial fraction of soil moisture (θg/θd, dimensionless). Where θg is the soil water content based on mass expressed in the gravimetric soil moisture content, θd is the soil moisture. Uppercase letters (V, U, M and G) represent the best calibrated simulations in the wet season and dry season of 2014, analyzed by means of Taylor diagrams.

Parameter		Wet Season Simulation V	Dry Season Simulation U	Wet Season Simulation M	Dry Season Simulation G	Source
	Initial	Energy Flux Calibrated		CO2 Flux Calibrated
Specific leaf area (sla)	13.0	14.5	23.5	14.5	23.5	[61,62]
Typical dimension of leaves (du)	0.072	0.056	0.032	0.056	0.032	[63]
Typical dimension of stems (ds)	0.1	0.05	0.05	0.05	0.05	[63]
Leaf width (w)	0.1	0.06	0.03	0.06	0.03	[63]
Coefficient of stomatal conductance (m)	10.0	8.0	5.0	8.0	5.0	[4,61,62]
Maximum capacity of the Rubisco enzyme (Vmax)	75 × 10−6	90 × 10−6	90 × 10−6	90 × 10−6	60 × 10−6	[62,64]
Initial fraction of soil moisture (θg/θd)	0.36	0.165	0.075	0.225	0.165	[38,39]

. Moreover, these parameters were selected not only because of their importance and sensitivity in driving the main components of the energy balance and CO₂ fluxes but also due to the high uncertainty related to their in situ measurements.

The model performance was evaluated and analyzed by means of Taylor diagrams [57] using the Pearson correlation coefficient (r), standard deviations (SD) between observations (x-axis) and simulated data (y-axis) and root mean square error (RMSE). Thus, SD > 1 indicates that the model overestimated the variables, and SD < 1 indicates that it underestimated them. In a Taylor diagram, radial distance represents the ratio of simulated to observed standard deviations, the azimuthal angle represents simulated-observed data correlation, and the distance between observed and simulated data points corresponds to the RMSE. These diagrams are useful because they provide summarized information on the relative performance of an ensemble of simulations. Furthermore, the model performance was assessed through the following statistical measures: median absolute error (MAE) and Willmott’s index of agreement (d). For plotting Taylor diagrams, we used the taylor.diagram function of the R package plotrix [58]. The daily means and totals of the fluxes variables were bootstrapped over seasonal intervals for the estimation of random variance (±95% of confidence interval—CI) about the mean according to the methodology presented in the literature [59]. Statistically significant differences (p < 0.05) in the mean seasonal value for a given meteorological variable or CO₂ flux components were determined by the degree of overlap in the 95% bootstrapped CI [60].

3. Results and Discussion

3.1. Calibration Test

For a more comprehensive verification of the simulations we elaborated Taylor diagrams, which can display the correlation, standard deviations and root mean square error between observations and multiple simulations in a single diagram. Figure 1 and Figure 2 shows the results of 44 simulations of site-scale energy balance components, as well as CO₂ fluxes for the wet season (blue letters) and dry season (red letters) of 2014, compared to observations (open circles).

Initial calibrations of morphological and physiological parameters were crucial to obtain a more realistic representation of the energy fluxes in the Caatinga biome in both wet and dry seasons. Main calibrated parameters were: specific leaf area (sla), typical dimension of leaves (d_u), leaf width (w), coefficient of stomatal conductance (m), maximum capacity of the Rubisco enzyme (V_max), and the initial fraction of soil moisture (θg/θd) (Table 2).

In the Taylor diagram we can observe that simulated LE and H are most sensitive to soil moisture fraction and the coefficient of stomatal conductance, mainly during the dry season (Figure 1, Table 2). For θg/θd = 0.165, the model provides a better fit between simulated LE and H and observed values in the wet season, while setting θg/θd = 0.075 results in a better fit during the dry season (Table 2). On the other hand, the NEE simulated by the model is most sensitive to the V_max parameter (Figure 2; Table 2). When using the lowest V_max value (40 μmol CO₂ m⁻² s⁻¹) the SITE model underestimated NEE, while when using the highest V_max value (120 μmol CO₂ m⁻² s⁻¹) it overestimated NEE. The most appropriated value for V_max was 90 μmol CO₂ m⁻² s⁻¹ in the wet season and 60 μmol CO₂ m⁻² s⁻¹ in the dry season (Table 2).

After analyzing the calibrated simulations through the Taylor’s diagram, we compared the ones with the higher correlation coefficient, lower RMSE and standard deviation for each season (Table 2) with non-calibrated simulations and EC observations. Table 2 shows the calibrated parameters and the best calibrated simulations in the wet season and dry season (represented by uppercase letters V, U, M and G), as analyzed by means of Taylor diagrams. For the analyses of validation, the Pearson correlation coefficient (r), median absolute error (MAE), root mean square error (RMSE) and Willmott’s index of agreement (d) was included between the modeled data with calibration by SITE model and the observation data for the years 2014 and 2015 (

Table 3. Results of the Pearson correlation coefficient (r), median absolute error (MAE), root mean square error (RMSE) and Willmott’s index of agreement (d) between eddy covariance observations of energy and CO₂ fluxes and simulations with the parameters set after calibration for the years 2014 and 2015.

Statistic
	2014				2015
	r	MAE	RMSE	d	r	MAE	RMSE	d
Energy flux
Rn (W m−2)	0.98	43.84	80.46	0.94	0.96	31.70	50.56	0.98
H (W m−2)	0.85	50.52	68.83	0.89	0.91	46.94	71.97	0.89
LE (W m−2)	0.69	29.97	65.99	0.72	0.71	24.25	53.48	0.74
G (W m−2)	0.90	9.90	13.66	0.92	0.90	11.07	16.23	0.91
CO2 flux
GPP (g C m−2 h−1)	0.82	1.24	1.53	0.86	0.91	1.38	2.05	0.79
NEE (kg C m−2 h−1)	0.84	1.99	2.25	0.83	0.81	1.74	2.00	0.80

).

3.2. Daily Variations and Seasonal Dynamics of Simulated Energy Fluxes

Figure 3 shows the daily cycle of observed energy fluxes (solid black line), non-calibrated simulated data (dashed blue line) and calibrated simulated data (dashed red line) for Rn, H, LE and G during the wet season (February to May 2014) and Figure 4 shows the daily cycle of energy fluxes during the dry season (August to October 2014). Regarding observed values one can note that the components behave according to the daytime behavior of solar radiation, with maximum values occurring between 11h and 12h (local time). During the wet season, peak H and LE values were of 230 W m⁻² and 218 W m⁻², which corresponded to approximately 40% of Rn (564 W m⁻²), respectively (Figure 3). During the dry season, however, most of Rn (peak value of 582 W m⁻²) is clearly converted to H (peak value of 404 W m⁻²). The peak LE values are very small (20 W m⁻²) and the values of G and LE are approximately equivalent (Figure 4).

The non-calibrated simulation severely underestimated daily H, as much as 40% in the wet season and 66% in the dry season (Figure 3B and Figure 4B; respectively). Additionally, daily LE were overestimated by approximately 45% during the wet season and by up to 8.5 times in the dry season (Figure 3C and Figure 4C; respectively). After the calibration of the biophysical parameters (Table 2), the model satisfactorily described the daily variations of Rn, H, LE and G (r > 0.8, d > 0.7; Table 3). In general, the adjusted Rn and G were satisfactory, with r > 0.9 and d > 0.9 (Table 3). However, LE was not satisfactorily simulated (r< 0.7) in the dry season (Figure 1 and Figure 4). The main weakness of the model was in relation to the simulation of peak hours (approximately 12:00 local time), where H was underestimated in the wet season (15%) and the dry season (16%) (Figure 3 and Figure 4). On the other hand, simulated LE presented a good adjustment (r > 0.8) with observed data in the wet season, but overestimated it by approximately 35% in the dry season (Figure 4). These results indicate a certain difficulty by the model to represent variations of energy fluxes at peak times, even after exhaustive calibration tests.

The differences found between observed and simulated LE and H values during the day may be due to the fact that the model does not consider the daily variations of some parameters inherent to the ecosystem, such as the fraction of soil moisture. In the Caatinga biome, soil moisture decreases during the day, as reported by others [39,65,66]. This may explain the slight overestimation of LE after 12h00 in the wet season (Figure 3C). Similar difficulties when adjusting the hourly data variability was reported in flux simulations using the SITE model in a tropical semi-deciduous forest in the southern Amazon Basin [53] and also using other models, such as the Noah-MP, over a forest site in the Amazon [67] and integrated biosphere simulator (IBIS) for a Brazilian semiarid region [68].

The mean monthly values of the energy fluxes are presented in Figure 5and Figure 6.

Table 4. Mean value of observed and simulated (values in parenthesis) energy balance components data: net radiation (R_n), sensible heat flux (H), latent heat flux (LE), soil heat flux (G) and CO₂ fluxes: gross primary production (GPP) and net ecosystem CO₂ exchange (NEE) during the wet season and dry season for the years 2014 and 2015.

Variable	2014		2015
	Wet	Dry	Wet	Dry
Energy flux
Rn (W m−2)	164.6(172.5)	168.7(167.4)	162.6(173.5)	174.6(170.9)
H (W m−2)	59.7(58.8)	113.7(91.9)	64.6(67.8)	120.0(94.4)
LE (W m−2)	71.5(81.8)	5.2(18.6)	48.1(56.9)	4.4(17.5)
G (W m−2)	−0.8(0.4)	5.4(4.9)	2.3(2.8)	6.5(4.9)
CO2 flux
GPP (g C m−2 h−1)	0.26(0.29)	0.11(0.17)	0.20(0.25)	0.10(0.12)
NEE (kg C ha−1 h−1)	−0.66(−0.64)	−0.25(−0.34)	−0.69(−0.63)	−0.49(−0.40)

Bold values are significantly different between observed and simulated values at the 0.05 level.

shows the mean seasonal values of observed (and simulated) for the years 2014 and 2015 of the components of the energy balance. Mean seasonal Rn values ranged from 164.6 W m⁻² (172.5 W m⁻²) and 162.6 W m⁻² (173.5 W m⁻²) in the wet season of 2014 and 2015 to 168.7 W m⁻² (167.4 W m⁻²) and 174.6 W m⁻² (170.9 W m⁻²) in the dry season of 2014 and 2015, in accordance with the annual variation of solar radiation.

The mean simulated H fluxes in the wet season (58.8 W m⁻² in 2014 and 67.8 W m⁻² in 2015) were similar to observed values (59.7 W m⁻² in 2014 and 64.6 W m⁻² in 2015) (Table 4). There was a statistically significant difference between the means of H fluxes during the dry period, in which the simulated mean value was 91.9 W m⁻² in 2014 and 94.4 W m⁻² in 2015 while the observed value was 113.7 W m⁻² and 120.0 W m⁻², respectively for 2014 and 2015 (Table 4). Minimum (maximum) H (LE) values occurred during the wet season, while maximum (minimum) values occurred in the dry season. This reduction in LE and increased in H during the dry period is attributed to lower water availability in the soil due to the absence of rainfall over the Caatinga biome. In addition, we verified a statistically significant difference between simulated and observed LE mean values in the dry period (Table 4). Table 3 shows that the error statistics of H and LE between eddy covariance observations and simulations with the parameters set after calibration for the years 2014 and 2015 (values in parenthesis) presented a MAE value of 50.52(46.94), RMSE 68.83(71.97), r 0.85(0.91) and d 0.89(0.89) for H, and MAE value of 29.97(24.25), RMSE 65.99(53.48), r 0.69(0.74) and d 0.72(0.74) for LE. In general, mean seasonal G values were low, less than 10.0 W m⁻² and with a monthly variability similar to that featured by H. The highest G values occurred during the dry season, period in which leaf senescence occurs in the Caatinga and the soil is more exposed to radiation. The analysis of G fluxes simulated by the model showed no statistically significant difference between simulated and observed values (Table 4).

Both the underestimation and overestimation of H and LE retrieved using the non-calibrated model (dashed blue lines, Figure 5) were greatly reduced with the adjustment of the SITE model according to specific wet season and dry season parameters (Figure 5 and Figure 6). The model was run twice, and the specific forecasts of each season were merged. These discrepancies in the non-calibrated model happened mainly because under conditions of low water availability a higher portion of the net solar radiation is converted to H. Our results were consistent with observational analysis in a tropical semi-deciduous forest in the southern Amazon Basin performed according to the literature [53]. These authors reported that LE was overestimated by 40% and H was underestimated by 66%. The discrepancy between observed and simulated H and LE values was observed in other studies using the SITE model in the Amazon rainforest [35] and also using other models in other ecosystems, such as temperate grasslands [69].

As discussed in the literature [70], throughout the annual cycle, LE behaves as a function of annual mean rainfall due to water limitations observed in tropical savanna ecosystems during the dry season. In addition, the energy flux partition is directly associated to vegetation characteristics and land use changes [16,33]. The occurrence of higher H values in arid and semiarid regions is a consequence of the reduction in water availability caused by low rainfall in these regions. Some studies show that in the Caatinga most of the available energy is converted into sensible heat flux in the dry season [33,34]. A similar trend of decreasing LE during the dry season was reported in the literature [34], which studied energy balance components in a preserved Caatinga region. These mechanisms of LE variability were satisfactorily simulated in the present study.

During the months with low water availability, the non-calibrated model greatly underestimates H. The default input parameters probably lead to the simulation of an ecosystem facing a higher water stress in relation to the Caatinga observational data. Thus, the use of specific parameters (e.g., reduced initial fraction of soil moisture—θg/θd) during the calibration of the model was necessary to improve its performance in simulating LE. The same issue when adjusting variations in H and LE values during the dry season was reported for simulations of energy fluxes using the SITE model in a semi-deciduous tropical forest in the Amazon [53].

The seasonality of rainfall and light greatly affect leaf morphology and physiology in the Caatinga species. During the dry season, the Caatinga vegetation, which consists primarily of deciduous and semi-deciduous species, reaches its minimum level of physiological activity, which in turn is barely sufficient for the maintenance of leaves. Indeed, leaves become shorter, narrower and smaller, leading to a significant decrease in the specific leaf area [61,63]. Under these conditions, transpiration rates in the Caatinga are very close zero [62,71]. The direct effect of these characteristics is the reduction of LE, with Rn being mostly converted into H.

3.3. Daily Variations and Seasonal Dynamics of Simulated CO₂ Fluxes

Figure 7 shows the daily variations in observed and simulated GPP and NEE, considering both non-calibrated simulations and the simulations with the best calibration for the wet and dry seasons. The observed GPP and NEE values presented a pronounced daily cycle, with larger amplitude during the wet season and smaller amplitude during the dry season (Figure 7), indicating that carbon uptake was more intense in the wet season due to the higher water availability. The observed GPP and NEE values sharply increased after sunrise, reaching maximum values between 09h and 11h and declining after noon, reaching its lowest values (near zero) at the end of the afternoon.

The non-calibrated simulations greatly overestimated the mean daily cycle values during both seasons (Figure 7). This illustrates the effect of neglecting the contribution of the physiological and morphological parameters of plants in this environment when predicting CO₂ fluxes. The main parameter influencing NEE and GPP was V_max. Using a value of 90 μmol CO₂ m⁻² s⁻¹ for the wet season and 60 μmol CO₂ m⁻² s⁻¹ for the dry season resulted in a fair agreement with observed values and a reasonable consistence with values obtained in previous ecophysiological studies involving species of the Caatinga [4,62,64].

In comparison with data from others research, we can observe that most of the modeling studies neglect the phenological variation of the ecosystems, usually considering a single value for V_max. For example, the values used in the Biosphere Energy Transfer and Hydrology model—BETHY [72] and in the global-scale mechanistic model HYBRID [73], 65 and 51 μmol CO₂ m⁻² s⁻¹, respectively, are similar to the dry season V_max used in the present study. In contrast, the values of V_max defined for the joint UK land environment simulator model—JULES [74,75] and community land model—CLM [76] are lower (48 and 31 μmol CO₂ m⁻² s⁻¹, respectively) if compared to V_max values used in other models. As reported in the literature [64] after calibrating V_max (approximately 90 μmol CO₂ m⁻² s⁻¹) values for the Caatinga in dynamic global vegetation models (DGVM), simulation results reached 72% of total observed GPP.

Simulations of GPP and NEE were also sensitive to changes in the coefficient of stomatal conductance, specific leaf area, typical dimension of leaves and leaf width (Table 2). The study [77] used a dimensional model for forest canopy radiation absorption, photosynthesis, and transpiration (MAESTRA) to estimate GPP in the Amazon Forest, showing that modeled GPP was sensitive to changes in total canopy leaf area and V_max, and that soil moisture, in addition to vapor pressure, controlled canopy CO₂ fluxes during drought periods.

The simulation also accurately retrieved the time at which the daily NEE peak occurs, between 09 h and 11 h. As presented in the previous study [35] using the SITE model, obtained a correlation coefficient of 0.88 between simulated and observed hourly CO₂ fluxes in a primary tropical evergreen forest, while in the literature [78] using two versions of the SSiB model in Amazon region, obtained correlation coefficients of 0.73 and 0.79 between observed and simulated hourly CO₂ fluxes.

The NEE is the difference between soil heterotrophic respiration and net primary production, thus negative NEE values indicate assimilation of carbon by the ecosystem. Furthermore, net primary production is GPP (gross photosynthesis) minus autotrophic respiration. NEE is highly dependent on stomatal opening, since stomatal conductance is commonly the main determining factor in CO₂ fixation [79,80]. In addition, biochemical models of photosynthesis (P_N) show that P_N represents the minimum value of two limiting factors: maximum carboxylation rates of Rubisco (V_max) and electron transport rate (J_max) [81]. Thus, the decline in GPP and NEE after 11h can be explained by the influence of environmental factors on the stomatal opening and biochemical parameters of photosynthesis such as kinetics of the Rubisco enzyme [4,80].

One noteworthy exception in the accuracy of the model refers to the nighttime period (Figure 7C,D). Values of CO₂ fluxes during nighttime were overestimated by 48% (wet season) and 56% (dry season), leading to higher RMSE (2.25) and MAE (1.99) of NEE relatively to GPP, which presented RMSE of 1.53 and MAE was 1.24 (Table 3). Have been previously reported [9] that the sensitivity of soil-surface CO₂ fluxes to volumetric water content above 10 cm is probably closely related to a high percentage of root and microbial biomasses in the 0–10 cm profile of tallgrass prairie. Therefore, the overestimation of simulated night CO₂ fluxes can be explained by the combination of root metabolism, root exudation, and, consequently, microbial activity in the rhizosphere.

Figure 8 and Figure 9 shows of monthly mean values of CO2 fluxes in 2014 and 2015, respectively, evidencing the existence of a clear seasonal variability. Table 4 shows means seasonal values of observed (and simulated) for the years 2014 and 2015 of GPP and NEE. The GPP and NEE increased at the onset of the wet season and reaching peak values in April 2014 and 2015, which GPP declined until reaching values lower than 0.2 g C m⁻² h⁻¹ in the dry season. Mean seasonal GPP in 2014 varied from 0.11 g C m⁻² h⁻¹ (0.17 g C m⁻² h⁻¹) (dry season) to 0.26 g C m⁻² h⁻¹ (0.29 g C m⁻² h⁻¹) in the wet season (Table 4). In 2015, GPP values presented a similar trend, ranging from 0.10 g C m⁻² h⁻¹ (0.12 g C m⁻² h⁻¹) (dry season) to 0.20 g C m⁻² h⁻¹ (0.25 g C m⁻² h⁻¹) in the wet season (Table 4). The simulation was able to capture the GPP and NEE seasonal variability. This result may indicate an efficiency of the SITE model in representing of gross primary production (GPP) and net ecosystem CO₂ exchange (NEE) dynamics in the Caatinga environment.

In the wet season, climate, canopy, and soil water conditions were optimal for the Caatinga vegetation. Thus, the GPP and NEE peaked and surface energy exchange was driven mainly by LE, while H was low. Photosynthetic activity decreased at the onset of the dry season because the soil was dry, leaf senescence occurred and H dominated surface energy exchange. In other words, water stress modulates the variations in canopy CO₂ fluxes throughout the seasons. The increase in rainfall rates in the wet season favors CO₂ uptake resulting in increased GPP and NEE in the Caatinga. In contrast, a reduction in soil moisture content induces stomatal closure and might affect some photosynthetic traits such as the Rubisco kinetics—V_max [82].

Furthermore, the effects of drought in the GPP and NEE could not be captured by the SITE model unless the calibration for the dry season was carried out. These results demonstrate that net CO₂ fluxes are very sensitive to the physiological processes that control surface energy exchange. Moreover, differences between seasonal observational and modeled carbon pools highlight the importance of phenology as an essential tool for understanding productivity in Caatinga. Others researches also report the importance of phenology in the dynamic of global vegetation models [83,84].

4. Conclusions

This study evaluated the performance of the SITE model when incorporated with parameters and input data consistent with in situ observations from the Caatinga biome, which is a seasonally dry tropical forest in the Brazilian semiarid region. In general, the SITE model exhibited reasonable applicability to simulate variations in CO₂ and energy fluxes. We believe that the SITE model could be used to simulate a satisfactory vegetation response if we take into consideration the remarkable phenological seasonality of the Caatinga. The LE flux was the output variable with the most unsatisfactory adjustment, overestimating observed values in the dry season. Furthermore, our ecophysiological approach offers the possibility to explore morphoanatomical and physiological mechanisms determinants of GPP and NEE in the Caatinga biome. Evaluating and improving the representation of the vegetation structure, dynamics, energy and carbon cycle of the Caatinga in vegetation models will help further develop our understanding on the impacts of land-use changes on regional and global carbon cycle.