Estimation of aboveground net primary productivity in secondary tropical dry forests using the Carnegie–Ames–Stanford approach (CASA) model

The study was conducted at Santa Rosa National Park (SRNP; 10°50'N, 85°37'W), Costa Rica (figure 1). The SRNP receives 1391 mm of annual rainfall and has a mean annual temperature of 25 °C (Kalacska et al 2004). The vegetation is drought deciduous, with a mixture of pastures and secondary patches in various stages of regeneration (Kalacska et al 2004). The canopy includes young forests with 80%–100% of woody plants being deciduous during the 6 month dry season (rainfall <100 mm: December–May), to semi-evergreen forests with a 30%–50% deciduous vegetation in older stages of succession (Arroyo-Mora et al 2005). Field data collected from previous studies in SRNP were used (Calvo-Alvarado et al 2012, Hilje et al 2015). Nine plots of 0.1 ha were sampled during 2007–2010 following a standard methodology (Kalacska et al 2004). These plots were established in early, intermediate, and late successional stages of 21, 32, and 50+ years of age respectively, with three plots in each category. The early stage of regeneration is composed of shrubs and small trees, with open areas and a single stratum of tree crowns. The intermediate stage is composed of deciduous trees and lianas, and it has two vertical strata. The late successional stage has two strata, and is formed by a dominant canopy layer and regeneration of shade tolerant species with reduced light penetration (Kalacska et al 2004, Arroyo-Mora et al 2005).

Zoom In Zoom Out Reset image size

Ground measured ANPP (ANPP_mea) between 2007 and 2010 was derived from litterfall production (Calvo-Alvarado et al 2012) and tree biomass increments (Hilje et al 2015). Biomass for each tree was calculated using a pan-tropical allometric equation that is a function of diameter at breast height, species-specific wood density, and an environmental variable denoting the dependence of bioclimatic changes (e.g., rainfall and temperature) on tree biomass (Chave et al 2014). Biomass increment from tree growth and litterfall production are reported in units of carbon by using a factor of 50% (Hughes et al 1999). Details about how ANPP_mea was estimated can be found in supplementary information (section 1.1).

2.3.1. CASA overview

The CASA model assumes that net primary production is proportional to APAR (Bloom et al 1985), which enables the calculation of ANPP at large scale using LUE. Based on Monteith (1972), the CASA model uses LUE (ε) as proxy, treating ANPP_CASA as function of APAR:

$$\begin{eqnarray}&&{\rm{ANPP}}\;=\;\varepsilon \cdot {\rm{APAR}},\end{eqnarray} \tag{ 1 }$$

where APAR is calculated by incident photosynthetically active radiation (iPAR) at canopy level and fraction of photosynthetically active radiation (FPAR) absorbed by vegetation canopy:

$$\begin{eqnarray}&&{\rm{APAR}}\;={\rm{iPAR}}\cdot {\rm{FPAR}}\end{eqnarray} \tag{ 2 }$$

and ε is calculated by maximum LUE (ε_{max_npp}, or maximum conversion efficiency) limited by two scalars denoting effects from water and temperature stress:

$$\begin{eqnarray}&&\varepsilon \;=\;{\varepsilon }_{\mathrm{max}\_\mathrm{npp}}\cdot {{W}}_{{\rm{scalar}}}\cdot {{T}}_{{\rm{scalar}}}\end{eqnarray} \tag{ 3 }$$

ε_{max_npp} is usually a biome-specified variable, representing the maximum ability of a particular biome to convert absorbed radiation into dry matter. Although, some PEMs calculate ANPP by subtracting autotrophic respiration (Ra) from GPP, the CASA model estimates ANPP directly, and incorporates photosynthesis part used for autotrophic respiration in the ε_{max_npp} term (Reich et al 2006).

2.3.2. Estimating CASA parameters

The CASA parameters to be estimated in this study include iPAR, FPAR, temperature and water scalars, and ε_{max_npp} (figure 2). In order to derive ANPP_CASA in each successional stage, we generated a forest succession map for the study area (figure 2) (table 1). Time series data including MODIS products (MOD04/MYD04, MOD06/MYD06, and MOD13Q1) (table 1) and meteorological data were processed according to 16 d composite to calculate 16 d total ANPP_CASA for 2002–2013.

Zoom In Zoom Out Reset image size

Table 1.  Remote sensing datasets used in the current study.

Role in this paper	Parameter	Dataset	Spatial resolution (m)
Succession mapping	—	HyMap	15
iPAR estimation	Instantaneous solar zenith	MOD04/MYD04	10 000
	Ångström exponent	MOD04/MYD04	10 000
	Cloud top pressure	MOD06/MYD06	5000
	Cloud optical thickness	MOD06/MYD06	10 000
FPAR estimation	NDVI	MOD13Q1	250
	EVI	MOD13Q1	250

We took advantage of a wireless sensor network (WSN) at the SRNP to continuously measure transmitted and absorbed PAR and meteorological variables. A WSN is a collection of independent nodes, each one measuring micro-meteorological variables that transmits such information via wireless to a data aggregator that in turn sends the data via cellular or satellite link to a cyberinfrastructure remote site where it is processed using analytical techniques (Pastorello et al 2011). WSN data (PAR) was collected in SRNP (mainly intermediate successional stage) from 06 March 2013 through 01 February 2015 (see supplementary section 1.2 for details) and were used to test the reliability of MODIS data in SRNP.

Specifically the main information used in this study was: instantaneous above canopy iPAR, PAR reflected by canopy, and PAR transmitted through canopy.

• (a)  
  iPAR: is the amount of solar radiation in visible wavelength (0.4–0.7 μm) that can be absorbed by green canopy through photosynthesis processes (supplementary section 1.3). This study estimated regional 16 d integrated iPAR (iPAR_est) in SRNP using MODIS products based on the PARcalc method proposed by Van Laake and Sanchez-Azofeifa (2004, 2005) with some simplifications. One of the most important simplifications of the method is to ignore the presence of clouds during the dry season. The iPAR_est was evaluated by WSN measured iPAR (iPAR_mea) at the time span of 6 March 2013–1 February 2015.
• (b)  
  FPAR: is the fraction of PAR absorbed by green vegetation. Long time series of FPAR (FPAR_est) were estimated from vegetation indices such as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) (Prince and Goward 1995, Running et al 2000, Xiao et al 2004, Li et al 2007). We also examined relationships between measured FPAR (FPAR_mea) by WSNs, MODIS NDVI, and MODIS EVI using linear and logarithmic models from 6 March 2013 to 16 November 2014. The best fitting model was selected to produce long time series of FPAR_est (see supplementary section 1.4).
• (c)  
  Maximum LUE (ε_{max_npp}). The biome specific ε_{max_npp} in the literature is derived by minimizing difference between ANPP_CASA and ANPP_mea (Field et al 1995, Ruimy et al 1999). Nevertheless, since TDFs in different successional stages present distinct species compositions (Kalacska et al 2004), we assigned a specific ε_{max_npp} for early, intermediate, and late successional stages, by calculating the ratio of ANPP_mea and estimated APAR $({{\rm{iPAR}}}_{{\rm{est}}}\times {{\rm{FPAR}}}_{{\rm{est}}})$ using environmental scalars in the permanent plots based on equations (1)–(3), as follows:

  $$\begin{eqnarray}&&{\varepsilon }_{\mathrm{max}\_\mathrm{npp}}=\displaystyle \frac{{{\rm{ANPP}}}_{{\rm{mea}}}}{{{\rm{iPAR}}}_{{\rm{est}}}\cdot {{\rm{FPAR}}}_{{\rm{est}}}\cdot {{W}}_{{\rm{scalar}}}\cdot {{T}}_{{\rm{scalar}}}}.\end{eqnarray} \tag{ 4 }$$
• (d)  
  Temperature and water scalars (supplementary section 1.5). Temperature and water scalars limit conversion of solar radiation to ANPP in green vegetation. Temperature scalar, T_scalar(x, t), explains two patterns of plant acclimation to temperature: under extreme temperature and seasonal temperature swing. Water scalar, W_scalar(x, t), is a function of estimated and potential evapotranspiration, expressing water deficit from 0.0 to 1.0, where 1 is water saturation. Meteorological data (temperature and precipitation) was obtained from the meteorological station at SRNP.
• (e)  
  Successional map. We produced a map of the successional stages in order to estimate ANPP_CASA for early, intermediate, late successional stages in the whole SRNP area. Identification of successional stages employed a hyperspectral sub-pixel mapping technique called multiple criteria spectral mixture analysis (MCSMA) (Cao et al 2015). This technique was applied to a Hyperspectral MAPper (HyMap) image (Cocks et al 1998) acquired on March 2005. MCSMA was first applied by Cao et al (2015) to map secondary TDF succession at SRNP, and it proved to be very efficient in dealing with spectral variability in TDFs (supplementary section 1.6).

Across successional stages, the intermediate stand had the greatest ANPP_mea with an average biomass increment of 6.4 ± 2.5 Mg C ha⁻¹ yr⁻¹, and mean litterfall of 3.1 ± 1.0 Mg C ha⁻¹ yr⁻¹. Biomass increments and litterfall in late succession averaged 10.5 ± 1.2 Mg C ha⁻¹ yr⁻¹ and 5.9 ± 1.2 Mg C ha⁻¹ yr⁻¹, respectively. The early stage had a much lower biomass increment (5.2 ± 4.2 Mg C ha⁻¹ yr⁻¹) and litterfall (1.9 ± 1.5 Mg C ha⁻¹ yr⁻¹) (table 2). On average ANPP_mea at SRNP was 7.1 Mg C ha⁻¹ yr⁻¹, with rates of carbon uptake of 3.6 Mg C ha⁻¹ yr⁻¹ in early stages, 9.5 Mg C ha⁻¹ yr⁻¹ in intermediate, and 8.2 Mg C ha⁻¹ yr⁻¹ in late succession.

3.2.1. Estimated CASA parameters

3.2.1.1. Temporal dynamics and evaluation of estimated iPAR

The iPAR_est agreed robustly with iPAR_mea both in dry and wet seasons over the time series period (6 March 2013–1 February 2015) (figure 3(a)). The linear regression model also shows a strong relationship between iPAR_est and iPAR_mea (R² = 0.83, y = 0.997x) (figure 3(b)). For each iPAR_est, absolute errors ranged 0.6% to 14.8% with average of 5.0%. The maximum individual error of 14.8% arises at the transition time from wet season to dry season at 19 December 2013, with a significant presence of brown downed leaves. Because of the constant presence of clouds, both 16 d total iPAR_est and 16 d total iPAR_mea in wet season were lower than in dry season (115 × 10³ versus 159 × 10³ KJ m⁻²).

Zoom In Zoom Out Reset image size

3.2.1.2. Relationship between MODIS vegetation indices and ground measured FPAR

Figure 4 shows the NDVI-FPAR_mea and EVI-FPAR_mea relationship using both a linear and a logarithmic regression model. The FPAR_mea increased from 0.5 in the middle of the dry season to 0.95 in wet season (figure 4(a)). Although the correlation between EVI and FPAR_mea could be considered high with a R² = 0.83 (logarithmic model), versus R² = 0.76 (linear model) (figure 4(c)); the NDVI-FPAR_mea relationship appears similar for both models (logarithmic model, R² = 0.91; linear model, R² = 0.90, figure 4(b)). Thus, posterior estimations of the FPAR_est for the modeling of ANPP_CASA were based on the NDVI-FPAR_mea relationship instead of the EVI-FPAR_mea model.

Zoom In Zoom Out Reset image size

3.2.1.3. Maximum LUE in different successional stages

Table 2 presents ANPP_mea and estimated APAR · W_scalar·T_scalar for the 2007–2010 period, and the ε_{max_npp} for each successional stage. Estimated ε_{max_npp} showed similar pattern with ANPP_mea. The highest ε_{max_npp} was presented in intermediate stages (0.81 g C KJ⁻¹), followed by late (0.66 g C KJ⁻¹) and early (0.31 g C KJ⁻¹) successions.

3.2.1.4. Seasonal dynamics of temperature and water scalars

Figure 5 presents boxplot of seasonal dynamics for the 16 d meteorological data and corresponding CASA scalar time series for complete year at SRNP. Each box was generated by summarizing corresponding meteorological data (12 July 2002 through 11 July 2014). Temperatures were stable across year at near the optimum value of 26.6 °C, with no significant differences between wet season and dry season (figure 5(b)). The temperature scalar also maintained high values greater than 0.9 (figure 5(d)). Precipitation, however, showed the two seasons characteristic in TDFs (figure 5(a)). The 16 d total precipitation was maximized in wet season, with almost no precipitation recorded in dry season. Precipitation presented its greatest inter-annual variation during the wet season, with a maximum standard deviation of 22.3 cm in late September (from 273th to 289th day of the year; not shown in figures). The water scalar had a similar pattern with precipitation. The land surface at SRNP became gradually water stressed as the dry season proceeded, despite that being relieved by occasional precipitation for instance on 03 April 2009 (6.7 cm) (figure 5(c)).

Zoom In Zoom Out Reset image size

3.2.2. Total and seasonal ANPP_CASA in different successional stages

We produced successional maps of early, intermediate, and late successional forests at SRNP by using the MCSMA to derive ANPP_CASA for each stage. Figure 6(a) compares area of TDFs in each successional stage at SRNP against their average ANPP_CASA (2002–2013). At SRNP, early successional stages comprised 56% of the total area, followed by intermediate successional stages with 37%. In contrast, early and intermediate successional stages comprised 32% and 59% of the total ANPP_CASA respectively. Both fractions of area and ANPP_CASA in late successional stages were small. Figure 6(b) illustrates the variation of the ANPP_CASA in the dry and wet seasons (December–April; May–November). Each box (e.g., ANPP_CASA of early successional stages dry season) was generated using data from 2002 to 2013. For all successional stages, ANPP_CASA in the dry season (early: 1.06 Mg C ha⁻¹ yr⁻¹; intermediate: 3.04 Mg C ha⁻¹ yr⁻¹; late: 2.67 Mg C ha⁻¹ yr⁻¹) was half of ANPP_CASA in the wet season (early: 2.17 Mg C ha⁻¹ yr⁻¹; intermediate: 5.86 Mg C ha⁻¹ yr⁻¹; late: 4.91 Mg C ha⁻¹ yr⁻¹). For each year, the ANPP_CASA of intermediate stages (8.90 Mg C ha⁻¹ yr⁻¹) was higher than ANPP_CASA of late successional stages (7.59 Mg C ha⁻¹ yr⁻¹), and 2.8 times higher than ANPP_CASA of early successional stages (3.22 Mg C ha⁻¹ yr⁻¹) (figure 6).

Zoom In Zoom Out Reset image size

Our estimates of ANPP_mea from ground data are similar to another study conducted in Costa Rica' TDFs. Waring et al (2015) reported rates of carbon gain from biomass increments and litterfall production 4.04 Mg C ha⁻¹ yr⁻¹ to 5.73 Mg C ha⁻¹ yr⁻¹ on stands of 15–65 years of age since land abandonment. For TDFs in other regions, Lugo and Murphy (1986) estimated ANPP_mea in Puerto Rico's Guanica forest as 3.45 Mg C ha⁻¹ yr⁻¹ by summing forest biomass increments and litterfall production. Martinez-Yrizar et al (1996) reported ANPP_mea of 3.06, 3.14, and 4.04 Mg C ha⁻¹ yr⁻¹ at three plots with decreasing elevations within a watershed at Chamela Biological Station, Mexico. Their measurements, however, included other two ANPP components leaf herbivory and understory production. When only biomass increments and litterfall production were considered, ANPP_mea for the three plots at Chamela Biological Station were 2.65, 2.74, and 3.56 Mg C ha⁻¹ yr⁻¹, respectively, which are comparable to our measurements. The variation on the upper limits of ANPP among existing studies and with our results could be explained by higher water availability at SRNP (mean annual precipitation; Guanica: 860 mm yr⁻¹; Chamela: 707 mm yr⁻¹; SRNP: 1390.8 mm yr⁻¹).

The relative lower correlation between EVI and FPAR_mea could be partially a result of the high sensitivity of MODIS EVI to the Sun-sensor geometry effect and canopy structure (Morton et al 2014). However, other studies have found better correlations between MODIS NDVI and LAI than using MODIS EVI in deciduous forests (Wang et al 2005) and in TDFs. Silveira et al (2007) for example, found that the best vegetation index for mapping vegetation classes in deciduous and semi-deciduous forests and Cerrado (Brazilian savannas) were obtained using the MODIS NDVI images than using MODIS EVI. This may be due to the fact that EVI tends to be more sensitive to NIR reflectance (Huete et al 1997) since it is more responsive to canopy structural variations, including LAI, canopy type, and canopy architecture (Huete et al 2002), whereas the NDVI is more chlorophyll sensitive (Huete et al 2002). This might be the reason why NDVI works better in dry forests compared to wet forests, since dry forests have a very heterogeneous forest canopy, with greater canopy openness and lower LAI even in the wet season (Arroyo-Mora et al 2005, Kalacska et al 2005b, Castillo-Núñez et al 2011).

Meteorological conditions dominate seasonal ANPP patterns at SRNP. Temperatures were stable at near optimum values across the year, making it an insignificant factor in ANPP_CASA estimation. Precipitation and FPAR exhibited opposite seasonal patterns to iPAR because of occurrences of rainfall and cloud cover. Precipitation greatly varied between dry and wet seasons. It is not surprising that abundant precipitation in the wet season promotes high ANPP_CASA because water availability is one of the main controls of leaf production and photosynthesis in TDFs (Jaramillo et al 2011). Our ANPP_CASA estimations showed TDFs sustaining photosynthesis even in the driest months at SRNP. This is probably related with changes in species composition across successional stages, with more than 80% of plants losing their leaves in the dry season, while only 30%–50% deciduous species in older stages of succession (Kalacska et al 2004, Arroyo-Mora et al 2005). Several plant species in TDFs have evolved different adaptive mechanisms with deep roots ensuring water supply for photosynthesis (Nepstad et al 1994). For example, woody vines, which are specially abundant on intermediate TDFs stages of succession, uptake more water than trees during water stress periods (Chen et al 2015) and as such tend to drop their leaves later in the dry season (Kalacska et al 2005a).

FPAR, which expresses the forest canopy structure and greenness (Arroyo-Mora et al 2005, Kalacska et al 2007), had a similar seasonal pattern with the water scalar in TDFs with high values in wet season and low values in dry season. This confirms the importance of rainfall in the leaf phenology of TDFs. The cyclical regimes of precipitation largely drive leaf flushing and falling events in secondary TDFs across different latitudes (Martha et al 2013). In our study, seasonal variability of ANPP_CASA from iPAR was dominated by FPAR and the water scalar, and values of ANPP_CASA during wet season were twice as much as in dry season. This differs from other studies in tropical environments that concluded that iPAR was the most influential climatic factor for primary productivity (Imoto et al 2010).

Assuming that SRNP experienced homogeneous meteorological conditions (iPAR_est, temperature, precipitation) across study plots, variations in ANPP_CASA originated from differences in FPAR_est and ε_{max_npp} across successional stages (see equations (1)–(3)), likely explained by differences in species composition and stem density. At tree level, ε_{max_npp} is a function of tree physiological processes and allometry, and in turn, a function of tree species. TDFs in late successional stages are dominated by shade-tolerant species with lower growth rates, while TDFs in early and intermediate stages have greater abundance of pioneer species that prefer full sunlight conditions with faster growth rates (Carvajal-Vanegas and Calvo-Alvarado 2013). As a result, early and intermediate stages present higher tree diameter increments (early: 1.6 mm tree⁻¹ yr⁻¹ versus intermediate: 2.2 mm tree⁻¹ yr⁻¹) than in late successional stages (1.2 mm tree⁻¹ yr⁻¹) (Carvajal-Vanegas and Calvo-Alvarado 2013). At regional level, ε_{max_npp} is further a function of stem density (Kalacska et al 2005b), promoting TDFs to reach their highest ANPP (mainly by litterfall production) in intermediate and late stages (table 2). This highlights the important role of intermediate and late successional stages in carbon sequestration, since these two successional stages have a greater ability to convert absorbed solar radiation into plant primary productivity (intermediate: 0.81 g C KJ⁻¹ versus late: 0.66 g C KJ⁻¹) than early successional stage (0.31 g C KJ⁻¹).

Furthermore, FPAR reflects canopy differences between different successional stages. Older successional stages with lower canopy openness and deciduousness have higher FPAR than younger successional stages (Arroyo-Mora et al 2005). For instance, TDFs in early succession were dominated by short trees, shrubs and grasses, which translated into higher canopy openness and lower greenness (Sanchez-Azofeifa et al 2009), and thus had lower FPAR_est compared to intermediate and late successional stages. It has also been reported that open canopies in TDFs are more vulnerable to wind and storm effects, losing their leaves faster and sooner in dry season than closer canopies (Jaramillo et al 2011, Calvo-Alvarado et al 2012). The successional effects from ε_{max_npp} and FPAR_est together explained that, early successional stages that currently represented more than half of the total area of SRNP accounted for only one third of the total ANPP_CASA, indicating that TDFs ANPP are more influenced by forest succession and species composition rather than forest area and forest extent.

The dominating role of precipitation in seasonal ANPP variation of TDFs highlights the need to collect more accurate and spatially explicit precipitation or soil moisture data in PEMs, despite that only one meteorological station was available at SRNP for our study. It is also important to consider the influence of structure and species composition of forest stands on ANPP in TDFs. Our model indicates that species composition explained the variation of ANPP. Other studies in TDFs in Costa Rica have found that these changes in species composition across successional stages may be explained not only by previous land use (e.g., stand age), but also by soil properties, including soil moisture (Becknell and Powers 2014). The water availability in our model relies mostly on rainfall, and we still lack a thorough understanding about the role of soil properties on ANPP. Future studies should explore the direct and indirect effects of soil on ANPP, via changes in species composition or by assessing the seasonal variation in soil moisture and its potential impact on rates of ANPP.