Modelling aboveground net primary production (ANPP) of an Atlantic mountain grassland based on time series approach

Because primary productivity is related both with the energy that sustains food webs and with species diversity, it is usually considered a key ecosystem property and a reliable indicator of available forage. In this work the aboveground net primary production (ANPP) of an Atlantic mountain grassland system was modelled in order to attempt producing short-term forecasts. Since grazing influences productivity, two treatment levels (grazing and exclusion) were experimentally applied in each of three field sites. Monthly ANPP data were then collected over three consecutive vegetative periods (2006-2008), thereby obtaining six time series (one per plot). Since no significant differences among sites (within treatments) were found, these six series were later reduced through averaging to only two series (one per treatment level). Two kinds of statistical models were then used to attempt monthly ANPP forecasting: exponential smoothing methods and ARIMA models. Both methodologies turned out to produce inadequate forecasts due to the presence of marked local features (innovative outliers) in our relatively short time-series data. Nonetheless, useful information for a more innovative shepherding management was revealed (e.g. the presence of within-year variation in ANPP, and differences between the grazing and exclusion treatments). Longer data series, which would require a more demanding effort in sampling investment, are likely necessary in order to obtain adequate forecasts using these time series methodologies.


Introduction
Plant productivity is, in essence, the dry matter produced by plants during a given time (Singh et al., 1975), where the underlying process is the fixation of atmospheric inorganic carbon into organic compounds by the photosynthesis. More specifically, aboveground net primary production (ANPP) is the fraction of gross productivity effectively fixed by aboveground plant organs once loses due to plant respiration have been taken into account.
ANPP is a key property of ecosystems because it is related with the energy that sustains the whole food web (McNaughton, 1985;McNaughton et al., 1989). Furthermore, some studies have suggested that the more energy available in the autotrophic level, the higher species diversity there is (De Angelis, 1980;Wright, 1983;Duffy et al., 2017), although this is a controversial topic and many others have asserted that there is no bivariate correlation between productivity and diversity Fox, 2013) or that this relationship strongly depends on the type of community (Mittelbach et al., 2001). In addition to being useful for community characterisation, ANPP is a reliable indicator of available forage and modulates the effect of grazing on several ecosystem functions (Díaz et al., 2007;Peco et al., 2017). For this reason, some studies have measured livestock carrying capacity based on productivity (Golluscio et al., 2015;Zhang et al., 2016), thereby providing a useful tool for management.
Productivity depends mainly on chemical soil properties (Fay et al., 2015;LeBauer and Treseder, 2008) and climate variables (Chapin III et al., 2011). Although annual precipitation is revealed as good predictor of broad-scale ANPP (La Pierre et al., 2016), at a local scale productivity is better explained through precipitation patterns and nutrient availability (La Pierre et al., 2016). Moreover, climate patterns which take into account growth cycles and phenology help explaining ANPP (La Pierre et al., 2011). Since ANPP influencing factors vary within years, ANPP is also thought to vary seasonally, except in tropical areas where seasonality may not be marked.
Some of the well-known consequences of the abandonment of livestock in grasslands include biomass accumulation (Patton et al., 2007), reduction of forage quality (Semmartin et al., 2004) and encroachment by woody plant species (Lasanta-Martínez et al., 2005;Gartzia et al., 2014). However, livestock grazing may affect grassland functions and ANPP via more subtle mechanisms. Grazing itself comprises not only the cutting and consumption of organic matter by grazers, but also trampling and fertilization by faeces deposition. For instance, low to moderate grazing intensity may enhance ANPP through grassland fertilization by animal urine and faeces (Frank and McNaughton, 1993;McNaughton et al., 1997), but it also may reduce ANPP via ecological mechanisms (non-specific disturbance) that allow the co-existence of highly productive species with less productive ones (Altesor et al., 2005;Tang et al., 2017). Moreover, the topsoil thermal regime and moisture conditions vary due to grazers activity (Gass and Binkley, 2011;Schrama et al., 2013), or because of livestock exclusion (Aalto et al., 2013;Odriozola et al., 2014). Therefore, ANPP fluctuates depending on the stocking density.
Despite the strong evidence suggesting variations in available forage due to changes in primary productivity (Fabricante et al., 2009), few studies have investigated withinyear productivity fluctuations. This is likely due to methodological limitations leading to lack of appropriate data, because primary productivity is ordinarily estimated as annual productivity rather than as monthly productivity (La Pierre et al., 2011Pierre et al., , 2016. To the best of our knowledge, previous attempts to model ANPP have not tackled this within-year variation by means of time-series modelling and, for this reason, they may have missed essential features in data structure useful for the prediction of plant biomass growth. To explore this issue, we conducted an experiment in the Aralar Natural Park, using monthly (not annual) time series data of ANPP, together with an appropriate set of statistical techniques such as time series data decomposition and forecasting. We aimed to measure and model change in primary productivity of Aralar grasslands under either grazing or grazing exclusion conditions, and to attempt the production of short-term forecasts.

Study area
The Atlantic grasslands of the Basque mountains provide important food resources to livestock, managed in short transhumance, from May (or June) to October (or November). Our experiment was developed in the semi-natural temperate grassland system of the Aralar Natural Park (11,000 ha), Basque Country, Northern Iberian Peninsula. The bedrock in the Aralar range is mainly calcareous (Gibbons and Moreno, 2002). The climate is oceanic temperate, with average (2006)(2007)(2008) annual precipitation equal to 1331 L m -2 and mean monthly temperature in the same interval equal to 7.0ºC. Due to this rainfall regime, the upper layers in the Aralar soils have become moderately acidic. Traditionally, the Aralar area has been used for livestock grazing and this activity has shaped the landscape, creating and maintaining large prairies. As a consequence, graminoid species diversity is high, with Festuca nigrescens Lam. and Agrostis capillaris L. being the dominant species (Odriozola et al., 2017). These grasslands are phytosociologically called Jasiono laevis-Danthonietum decumbentis (Loidi, 1983) and are considered a priority habitat (6230 subtype a) by the Habitat Directive (EC, 2013). Altogether, the habitat occupies more than 1980 ha within the Special Area of Conservation Aralar (BOPV, 2016).
In our experiment, monthly aboveground net primary production (ANPP) was defined as the weight (dry matter) of live biomass grown during one month (Singh et al., 1975). Monthly ANPP was measured in both fenced and grazed plots between January 2006 and December 2008, harvesting three to four 1 m 2 -quadrats per plot. Successive aboveground cuttings of standing crop of biomass were carried out at intervals of 25-30 days during the vegetative period (from the beginning of May, just after the snowmelt, to late October, before winter), so six monthly ANPP were measured per year and plot during the growing season. Since the average temperature in the non-vegetative season (between November and April) was consistently lower than 5ºC, the productivities of those months were insignificant and hence assumed to be zero. Thereby 36-month productivity values were recorded in each of the three fenced plots and an equal number were recorded for the grazed plots. Therefore, the whole production cycle was covered. Plant standing biomass was clipped at the ground level to ensure a uniform clipping height, and the percentage of live biomass was estimated visually. Samples were later dried at 65ºC for as much time as required (at least 48h), so as to avoid moisture-dependent weight fluctuations. Finally, mean live biomass (grams of dry matter per square metre) was obtained per plot before calculating monthly ANPP.
Since there are different methods to calculate ANPP, we chose the positive increment of live biomass recommended by Singh et al. (1975). The procedure of ANPP estimation was slightly different depending on the treatment. In the fenced plots, mean live biomass (LB) was calculated by harvesting randomly chosen 3-4 quadrats (1 x 1 m). The ANPP of a known month could be measured by subtracting the mean live biomass (LB fence ) of the previous month (t n ) from the live biomass of the current month (t n+1 ). In the grazed plots, 3-4 temporary cages (1.5 x 1.5 m) were used to avoid livestock eating herbage for one month (Frank et al., 2002). Standing live biomass was then clipped simultaneously inside the cages (LB cage ) and in the grazed plots (LB grazed ). Cages were then randomly moved to other points within the grazed plot allowing the growth of biomass during the next month (t n+1 ). In the case of grazed treatment, ANPP was measured by subtracting the previous month live biomass on the grazed plot (LB grazed in t n ) from the current live biomass inside the cages (LB cage in t n+1 ). That is: Monthly ANPP from t n to tn+1 in fenced treatment: Mean LB fence, n+1 -Mean LB fence, n Monthly ANPP from t n to tn+1 in grazed treatment: Mean LB cage, n+1 -Mean LB grazed, n

Data analysis
In order to model monthly ANPP between January 2006 and December 2008, several statistical techniques for time series analysis were conducted using R software (R core team, 2015). Since ANOVA analysis of this very same dataset provided no evidence of between-site differences (see Supplementary Material, Appendix 1), average monthly ANPP series according to treatment were computed. Hence, the original six data series were reduced to just two: one (average) time series for the fenced plots and another one for the grazed plots, each of them composed of 36 (average) values.
In the first of our analyses the characterisation of both time series was achieved through decomposition into three components: trend, regular component (seasonality) and irregular component (noise). Both additive and multiplicative decompositions (Cryer and Chan, 2008) were attempted by means of the function decompose() (package "stats" of base R), because the goal of this analysis was to explore which type of model better fitted the available data, which exhibit a marked feature: an unexpected increase in productivity for 2008.
Productivity for the following years could then be forecasted through either explicative regression models (using variables such as temperature or precipitation as predictors) or non-explicative models. Since the aim of this study was to predict future values based on past values, two non-explicative analyses were conducted. After careful consideration of the potentially useful available models and the fact that a priori we ignored which could be the best, we attempted modelling using both the so-called exponential smoothing method and seasonal ARIMA (autoregressive integrated moving average) models. Therefore, in the second set of analyses, one exponential smoothing method was used (Coghlan, 2014) in order to produce shortterm ANPP forecasting. As decomposition revealed both linear trend and seasonal components, a model of type A-A (Hyndman and Khandakar, 2008) was found to be appropriate. However, the increase of seasonality in 2008 suggested that the data could better fit a multiplicative seasonal A-M model (Hyndman and Khandakar, 2008). Consequently, exponential smoothing methods of A-A and A-M series were applied. Function HoltWinters() (package "stats" of base R) was used to produce one-step ahead forecasts called traditionally Holt Winters methods (Winters, 1960;Holt, 2004).
In the third analysis, a model-building strategy (Cryer and Chan, 2008) was developed in order to produce seasonal ARIMA model forecasts (Hyndman and Athanasopoulos, 2012) for each treatment level. First of all, model assumptions were checked (Cryer and Chan, 2008) by means of time series plots, autocorrelation plots (ACF) and partial autocorrelation plots (PACF), created through the function tsdisplay() of the package "forecast" (Hyndman, 2017). A logarithmic transformation was applied to meet the criterion of additivity and seasonal differencing was subsequently applied to achieve stationarity. Once both series were suitable for applying ARIMA models, a tentative model was chosen through the function auto.arima() of the same package. In the next step, function arimax() of the package "TSA" (Chan and Ripley, 2012), which uses maximum likelihood, confirmed the need to keep in the fitted models all the parameters suggested by autoarima(). Outlier values were then modelled to improve overall performance, using detectIO() and detectAO() of "TSA" (Box et al., 2008). Finally, when the models were validated (i.e. after checking for normality, null mean and independence contrasts of the residual values), future values were forecasted using forecast.Arima() of "forecast".
The R code including the original dataset and the three analyses are given in the Supplementary Material, Appendix S1.

Characterisation by time series decomposition
Decomposition revealed an increasing slope, a feature that is essentially due to the high productivity observed in 2008 (Fig. 2, trend component). The seasonal coefficients revealed a big peak of productivity in May and a second smaller peak in August in both series (Fig. 2, seasonal component). A third peak at the end of summer appeared, but only in the fenced series (Fig. 2, seasonal component). Thereby, the monthly ANPP series multiplicative decomposition showed the presence of changes in productivity, both within years (seasonality) and between years (positive linear trend) (Fig. 2). The results were similar when the additive decomposition model was tried (Fig. S2 in Appendix 2, Supplementary Material).

Figure 2. Multiplicative decomposition of monthly ANPP measured in Aralar field plots. Productivity (dry matter in g m -2 per month) is plotted in green (averaged ANPP under grazing conditions) and in dashed red (averaged ANPP measured in the fenced plots). Decomposition shows larger seasonal components under exclusion conditions. Trend components show high monthly productivities in 2008
, also more markedly in exclusion. All components are given in g m -2 per month.

Exponential smoothing forecasts
As revealed by visual inspection and residual errors, both additive (Fig. 3a) and multiplicative (Fig. 3b) Holt-Winters models did not fit well the available data. Nonetheless, model fit was more acceptable for the grazed-level averaged series than for the fenced-level averaged series. For this reason, forecasting based in these models cannot be trusted, as revealed by the wide 95% prediction intervals.

Figure 3a. Additive Holt Winters methodology forecasts. Forecasted productivity (dry matter in g m -2 per month) is between mean values of 2006 and 2008. Only two peaks are predicted at the end of both series. The sum of square errors (SSE) is equal to 45194 for the fenced series,
whereas equal to 7157 for the grazed series.

ARIMA models forecasts
No autoregressive or moving average processes were found with the exception of a first order seasonal autoregressive process. An especial effort was carried out to model outliers, which were included in the fitted models to satisfy residual normality requirements (Table 1). Thus, ANPP forecasts were produced through the equations specified in Table 1. ARIMA models, particularly in the case of the grazed-level series, fitted better the available data than Holt-Winters models and thereby forecasted more credible ANPP values (Fig. 4). Nonetheless, the prediction interval for the fenced-level series was huge ( Fig. 4; Fig. S3 in Appendix 2, Supplementary Material).

Discussion
Understanding both between-and within-year variations in ANPP is essential since ANPP is a key ecosystem property related with the energy sustaining the whole food web (McNaughton, 1985;McNaughton et al., 1989), as well as to biodiversity (De Angelis, 1980;Wright, 1983;Duffy et al., 2017). Any advances on the prediction of ANPP variations at annual or seasonal temporal scales would have implications for management. ANPP interacts with grazing pressure affecting multiple ecosystem functions (Díaz et al., 2007;Peco et al., 2017), thus any knowledge on future ANPP variations would allow a better adaptation of livestock carrying capacity based on productivity (Golluscio et al., 2015;Zhang et al., 2016), with consequent optimisation of related ecosystem functions. Because of the short series used, our study failed on accurately forecasting monthly ANPP variation. However, we found interesting features for management related to seasonal and annual variations on ANPP, and most importantly, future studies on forecasting ANPP will benefit from our detailed discussion on different statistical methodologies implemented, together with the R code, made freely available on the Appendix 1 of this work.
In this work, ANPP, which was measured straight away after fencing, was almost twice that measured under grazing conditions. Since, in broad terms, similar or even higher ANPP are often reported under grazing conditions because of greater renewal of leaf tissue and lower senescence rates, the greater ANPP found in our study in excluded vs. grazed plots requires additional explanation. We believe that this is a short-term effect, most likely due to the observed fact that the growth of the dominant grasses (Festuca nigrescens and Agrostis capillaris) was released when large grazers were prevented by fencing (Odriozola et al., 2017). Compared to forbs, grasses strongly determine ANPP, and grazing may cut down the response of grasses to increased temperature (Wang et al., 2012). Besides, Patton et al. (2007) described negative effects of grazing abandonment on ANPP, but only after long-term grazing exclusion, when litter build-up inhibited the growth of the plants, and this may also be the case for Aralar. Both explanations contribute to support our suggestion that the observed difference in ANPP in excluded vs. grazed plots is temporary.

ANPP data collection
Although productivity is defined as biomass produced in a given time interval (Singh et al., 1975), there are many ways to collect biomass data and many time intervals to produce time series. Grassland productivity has been traditionally measured by nondestructive methods such as plate meter and 3D quadrat with varying success (Redjadj et al., 2012), or more reliable estimates under annual calibrations as radiometer and point frame methods (Byrne et al., 2011). Grassland productivity has been measured by more accurate destructive methods as well. Within the latter, ANPP has been underestimated considering the final peak biomass as the total production (La Pierre et al., 2011;Craine et al., 2012;La Pierre et al., 2016), or overestimated by not differencing between live and dead biomass (Hu et al., 2007). Additionally, studies have usually investigated productivity analysing it as annual (La Pierre et al., 2011 and not allowing the study of within-year variation. In this study, we used a destructive methodology, measuring monthly ANPP and differencing death and live biomass. This methodology allowed us measuring ANPP accurately, as well as investigating its seasonal variations. Nevertheless, this was achieved by a very time-demanding procedure. The development of less laborious but sufficiently accurate techniques would facilitate the collection of equally detailed but longer time series than ours.

Exponential smoothing forecasts
The Holt-Winters models did not fit well, in general terms, with the data: for the within-year peaks, the high productivity corresponding to year 2008 and the null productivity of winter are data features that were not described appropriately by the fitted models. As a consequence, both the sum of square errors (SSE) and the prediction intervals were too large. For those reasons, the application of this method on these data series turned out to be not satisfactory. Regarding the results of the multiplicative Holt Winters method, the same conclusion applies, as both the SSE and prediction intervals were enormous.
Since Holt-Winters methods forecasted different trends according to whether they were multiplicative or additive, further discussion is required so as to shed light on these differences. First of all, it is important to highlight that the increasing trend is due to the extraordinary ANPP values corresponding to 2008, which have influenced positively all future predicted values. If these high values had been located in the middle of the series, the trend component would not have been as influencing as it is. However, in view of the positive relationship of ANPP with either precipitation (La Pierre et al., 2016) or temperature (Craine et al., 2012), the presence of a positive linear trend, as in the additive method, appears to be logical, because in the period 2006-2008 a positive trend in temperature is apparent in the meteorological data. Nevertheless, in the fenced series there are other influencing ecological drivers than climate, because the community is changing due to the lack of equalizing mechanisms that herbivores exert, and the field plots were not in a stable equilibrium. In this case, a bigger increment of ANPP was revealed, which was related with the sudden spread of highly productive dominant grasses (Odriozola et al., 2017), and this enhanced the positive trend. As observed in long-term grazer exclusion experiments (Patton et al., 2007), the initial extra productivity would most likely cease after a longer period of grazing abandonment, when litter accumulation inhibits plant growth. Anyway, despite being the slope modified, a positive linear trend as forecasted by the additive method could be adequate for several years. By contrast, there is no evidence (variance increasing with trend component) to support the exponential positive trend suggested in the multiplicative method.
Summarising, additive Holt-Winters forecasts appear to be more trustworthy than multiplicative ones under grazing conditions, although additional studies recording longer time series data would no doubt be required to clarify trends under exclusion conditions. The multiplicative model, though, captured seasonality better, because productivity was positive in winter in the additive forecasts.

ARIMA model forecasts
ARIMA modelling did not lead to useful forecasting in either case, mainly for two reasons. On the one hand, the introduction of too many step functions to model outliers was inevitable in order to achieve a satisfactory model fit and reduce the residual errors. On the other hand, despite all these corrections, which inevitably introduced an unwanted mathematical complexity, the accuracy of the model was inadequate, a fact that is evidenced by the huge prediction intervals (especially under fenced conditions). These two reasons are, in turn, related with the shortness of the dataset. However, data variability in the grazed series was better captured using ARIMA modelling than using Holt Winters methods. In short, forecasts were not accurate and, once again, the collection of longer data series seems to be inevitable in order to obtain useful results from ARIMA modelling.

Implications for management
It is impossible to assess which of the methods is more appropriate for ANPP modelling, since both attempts have been unsuccessful. Nonetheless, although the main aim of this work (producing accurate forecasts) was not adequately satisfied, our analytical work and the script provided may be useful to apply in other longer-term datasets, or in combination with monitoring tools such as the one proposed by Primi et al. (2016). Additionally, our study did show other interesting features for grassland management. Seasonality in ANPP was proved, as in other studies which have linked it with the precipitation regime (Swemmer et al., 2007;Fabricante et al., 2009) or with plant phenology (La Pierre et al., 2011). This seasonality must, therefore, be considered at the time of planning the livestock management in the park. This will help to optimize the grassland utilization since it allows designing a temporal planning of the grazing so that livestock take advantage of the maximum peaks of forage.
As in many European grasslands (Pe'er et al., 2014) the reduction of livestock has been great in the mountains of the Cantabrian-Atlantic region, and the risk of woody plant encroachment due to the accumulation of biomass has been revealed (Gartzia et al., 2014;Lasanta-Martínez et al., 2005). Moreover, the decreasing trend in livestock is expected to accentuate in the future (Rounsevell et al., 2006). These grasslands with high conservation value (EC, 2013) are deeply linked to the action of herbivores. Experimental grazing exclusion has demonstrated that the decoupling of grassland ANPP and the livestock pressure has many unwanted consequences for their conservation, apart from shrub encroachment. These changes include loss of plant diversity (Odriozola et al., 2017); reduction on forage quality as reduced digestibility and nutrient content, and enhanced fibre content (Odriozola et al., 2014); and the related reduction on microbial activity and retarded nutrient cycling (Aldezabal et al., 2015).
In contrast to the current situation, other simulations of recovery of shepherding could be hypothesized. In this case, features also supported by other studies, like seasonality of the ANPP (Fabricante et al., 2009;La Pierre et al., 2011) and the increasing trend of the production (Craine et al., 2012;La Pierre et al., 2016), can help to increase the carrying capacity of the grasslands of the Aralar Natural Park. Active shepherding based on rotation systems might be used to create short-term abandoned fields, where forage production is enhanced. Moreover, the sheep flocks could also be actively guided to make better use of within-year (seasonal) productivity peaks.

Conclusions
The aboveground net primary production (ANPP) of an Atlantic mountain grassland system was modelled in order to attempt producing short-term forecasts. Monthly ANPP data were collected over three consecutive vegetative periods under grazing and exclusion conditions. Despite using several analytical methodologies, all of our attempts to forecast monthly ANPP accurately failed. This was mainly due to the relative shortness of the time series (n = 36 data points) and to the presence of numerous outliers, which has compelled us to introduce an undesired mathematical complexity in the ARIMA models. Nonetheless, our work revealed within-year variation in monthly ANPP and other useful information for the management of Cantabrian-Atlantic grasslands, such as the increase in productivity during the immediate years following grazing abandonment. Longer data collection should be necessary in order to use any time series forecasting methodology, although this would obviously imply more demanding inversion efforts.