Response of Seasonal Vegetation Dynamics to Climatic Constraints in Northeastern Burundi

Pacifique Batungwanayo^1,2*, Marnik Vanclooster³, Alan F. Koropitan⁴

¹Departement of Mathematics, Institute of Applied Pedagogy, University of Burundi, Bujumbura, Burundi.
²Research Center for Mathematical and Physical Sciences, Doctoral School, University of Burundi, Bujumbura, Burundi.
³Earth and Life Institute, Université Catholique de Louvain, Louvain-la-Neuve, Belgium.
⁴Faculty of Fisheries and Marine Sciences, IPB University, Bogor, Indonesia.

DOI: 10.4236/gep.2020.89010                

Abstract

The climate crises in East Africa (EA), particularly in Burundi, have affected vegetation which, in turn, plays a key role in the climate system by modifying the terrestrial water and energy balance. Consequently, it is vital to understand vegetation dynamics and its response to current and projected climate conditions to support the design of climate resilient land management strategies. The objective of this study was to study the dynamics of vegetation cover over the Northeastern Burundi (NEB) in response to climatic constraints. The methodology used consisted of the interpretation of satellite images along with the analysis of data collected through rain-gauge stations. The data sets used include time series composite moderate resolution imaging spectroradiometer (MODIS) normalized difference vegetation index (NDVI) data collected between February 2000 and December 2017; long term (1986-2017) rainfall data acquired from two meteorological stations throughout the Northeastern provinces in Burundi and precipitation and mean temperature data (1986-2017) from Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) and ERA5 Daily aggregates for the study area. The study provides an assessment of the vegetation trends in NEB using the NDVI time series. In addition, regression analysis is applied to assess the relations between NDVI and precipitation, air temperature, potential evapotranspiration, as well as with growing season characteristics (onset, cessation, and length). Results show that the rate in vegetation productivity is persistently gradual between 2000 and 2011 despite fluctuations from the mean position, followed by a lower growth rate over the period 2011-2017. There has been trend variation in precipitation, neither the temperature was constant. The temperature over the region has increased while the precipitation has decreased. The onset of the growing season and air temperature also show a significant influence on seasonal vegetation dynamics in the region. Drought-induced plant stress observed at the onset of the rainy season was the most important contributor to the subsequent less greening of vegetation especially over the area near the northern lakes.

Keywords

Rainfall Variability, Onset, Dry Spell, NDVI, Remote Sensing, Northeastern Burundi

1. Introduction

There is more evidence that climate change and variability are undeniable phenomena that affect the socio-economic activities of many countries all over the world, including Burundi. It has received more attention at the global level and has serious consequences for food security and livelihoods, particularly in developing countries (UNFCCC, 2007; Manyeruke et al., 2013; Tirivangasi, 2018). In most parts of Sub-Saharan Africa (SSA) for instance, the main farming practice is rainfed agriculture (Beyer et al., 2016). About 97% of cultivated land in SSA is dominated by rain-fed agriculture (Calzadilla et al., 2013) and this farming practice is expected to remain the main source of basic food production for most rural communities (Cooper et al., 2008). As developing countries are highly dependent on agriculture, there is a growing concern that a shift in climate variability will further threaten the well-being and food security of already highly vulnerable rural households, leading therefore to serious challenges for development efforts (Abdulraheem et al., 2017). Despite the efforts made in advanced technologies supporting food security such as improved varieties, genetically modified organisms, and irrigation systems, climate exerts major pressure on agricultural productivity and hence food security. The effect of climate on agriculture is principally related to the variability of local climates rather than the global climatic regimes (Field et al., 2017). Hence, assessing the variability of climate properties at the local scale is pivotal for increasing food security. This implies that the study of the spatiotemporal variability of rainfall, temperature, and humidity at the earth’s surface at the local scale is an important prerequisite for designing sustainable agricultural and food development programs.

With increasing climate variability and warming trends in SSA, particularly in EA, there has been a shift towards extreme rainfall and an increase in mean seasonal temperature. Hence, these trends may have an impact on vegetation change and agricultural production, particularly in countries such as Burundi, where at least 90% of the population depends on rain-fed agriculture. There is however a need to better assess these possible impacts to mitigate to the effects of change and to design climate resilient food production programs.

In the last decades, Burundi went through periods of famine that raged and threatened thousands of people and still engrave in the memory of many Burundians (Ndayiragije et al., 2017). These include, for example, the famines in 1917, 1923, 1931, 1933, 1943 and 1958 throughout the country, 1989-1990 in the regions of Bujumbura Rural, Bururi, Gitega and Muramvya and 2000-2005 in the north-east of the country. Some of these famines were caused by adverse weather events (e.g. drought, excessive rainfall, and hail), whereas others resulted from calamities such as the invasion of predator insects (e.g. migrating grasshoppers and locusts) (Ndayiragije et al., 2017). Official recordings of weather-related perils in crop and livestock production are sparse in Burundi. Most prominent weather-related calamities were prolonged rainfall deficits and torrential rains. Future climate predictions show that precipitation for Burundi is projected to increase in the eastern and southern parts of the country as well as in the central plateau. While rainfall will increase during the rainy season, the months before the onset of the rainy season (August up to September) may become drier and longer (Liersch et al., 2014). This may result in postponed planting dates and ultimately may cause harvest losses. There is also a high probability that annual average air temperatures will steadily increase in Burundi over the century. Liersch et al. (2014) assert that the changes in precipitation patterns and quantity, as well as temperature may have significant implications for crop production.

The Northeastern part of Burundi is one of the most sensitive regions to changes in temperature, rainfall, and water resources variabilities in the country involving food shortage and insecurity. However, the effects of climate change on water availability and agriculture in the region have not been dealt with adequately, but are necessary to further understand the future impact of climate change on the region and to implement sustainable coping strategies. Only few studies related to an assessment of vegetation dynamics in response to climatic factors in order to mitigate their impacts have been done. However, most of the studies on how vegetation is responding to recent changes in climatic conditions are done at the regional scale over EA (Nicholson et al., 1990; Pelkey et al., 2000; Measho et al., 2019). There is no specific case studied for Burundi. It is needed to understand the vegetation phenological changes for planning purposes that can consequently reduce the loss of lives and food crisis over NEB. Therefore, a deep assessment of changes in precipitation and temperature and its impacts on vegetation dynamics is pivotal to mitigate climate change-related hazards and vulnerability. Such an assessment will help to implement climate-related policies.

In recent years, the use of remote sensing technology has increased globally as well as regionally to study vegetation dynamics and trends. The NDVI derived from satellite data has often been used to analyze the growing conditions of living green vegetation and reveal the response of vegetation dynamics to climate change (Yu et al., 2003; Measho et al., 2019). For EA for instance, Nicholson et al. (1990) assessed the relationships between NDVI and rainfall and found that spatial patterns of annually integrated NDVI closely reflect mean annual rainfall. Also for EA, Pelkey et al. (2000) observed an increase in greenness between 1982 and 1994. Considering the existent research situation, the behavior of the vegetation data series has to be well characterized in the NEB. Thus, the current study aims at filling the knowledge gap by assessing how vegetation in NEB is responding to recent changes and variability in climatic conditions. The specific objectives are to: 1) characterize the rainfall fluctuations on the both annual and seasonal time scale; 2) analyze the seasonal and yearly trends in NDVI time series; 3) extract phenological parameters, climatic parameters for short and long rainy seasons and 4) identify the NDVI-climate factors relationship on the both yearly and seasonal time scale.

2. Materials and Methods

2.1. Area of Study

Figure 1 shows the study area location in Burundi, an East African country ranging between longitudes 28.8˚E and 30.9˚E and latitudes 2.3˚S and 4.45˚S (Bidou et al., 1991).

Burundi, one of the world’s least developed countries, landlocked, located at 1200 km from the Indian Ocean and 2000 km from the Atlantic Ocean between

(a) (b) (c)

Figure 1. Area of study: (a) The locations of selected stations on the NEB region; (b) Mean monthly rainfall distribution from Kirundo station data; (c) Mean monthly rainfall distribution from Muyinga station data.

the DRC in the West, Rwanda in the North, and Tanzania in the East and the South, covers an area of 27,834 km². The country is ranked among the most densely populated regions of Africa and almost 90% of the population living in rural areas are involved mainly in agricultural activity (Niragira et al., 2015; Dinku et al., 2007). Although the country is in the equatorial climate regime, its high altitude leads it to a temperate climate that varies according to the topography. The NEB, consisting of the central plateaus and the northern depressions, experiences annual rainfall ranging from 1100 to 1500 mm. Covering about 13% of the country’s surface area, the study region provides on average, over its whole annual production, 30% of cereals, between 20% - 30% of legumes, 15% - 20% of bananas and plantains, and 10% of tubers and roots. This region has drawn particular attention in recent decades due to the fact that farmers experienced cropping losses, which led to hunger and immigration to neighboring provinces. Climate variability is cited as one of the reasons they could no longer farm as they used to (Ministry for Land Management and Environment, 2007).

The NEB experience two rainy seasons, season A running from Mid-September to December (SOND) known as the short rainy season and season B running from March to May (MAM) known as main rainy season. Both seasons are separated by two dry seasons, season C running from June to August (JJA) known as long dry season and a short dry season during January-February (JF). Season A represents on average one third of annual production while seasons B and C represent about 50% and 15% of annual production, respectively. The Inter-Tropical Convergence Zone (ITCZ) crosses the country twice a year and is one of the principal factors determining the rainfall seasonality (Nkunzimana et al., 2019).

2.2. Datasets

The study used multiple datasets from different sources to analyze the rainfall variability and its relationship with vegetation dynamics. The first set consists of observed data collected from the synoptic stations of Burundi belonging to the Geographical Institute of Burundi (IGEBU). Two meteorological stations have been considered over the historical period 1986-2017 and the concise information on their location is presented in Table 1. Since the African rain-gauge data have many spatial and temporal discontinuities over large sections of East Africa (Schreck III & Semazzi, 2004), the Google Earth Engine (GEE), an alternative platform for satellite data collection, has been used to get a second set of data. These include the CHIRPS data which incorporates 0.05˚ resolution satellite imagery with in-situ station data (Funk et al., 2015) and the fifth-generation ECMWF atmospheric reanalysis data of the global climate (ERA5). This latter dataset replaced its predecessor ERA-Interim reanalysis. The FLDAS dataset (McNally et al., 2017) designed to assist with food security assessments in data-sparse developing country settings that include evapotranspiration was also used in the study.

Table 1. Geographical description of two meteorological stations and length of climatic records considered in the study.

Furthermore, to assess the impact of rainfall variability on vegetation, the NDVI was generated from the MODIS/006/MCD43A4 surface reflectance composites for rainfall season on 16 days temporal resolution and 500 meters spatial resolution. The satellite imagery presents advantages of regularity in time and homogeneity in spatial coverage. All datasets were used to perform a comprehensive and detailed analysis in order to fill any gaps that station data may contain.

2.3. Methodology

The study aimed at investigating first the main features of rainfall over the study domain and, subsequently, the relationship between main variables that characterize the rainfall season (date of start of the short rainy season, date of the end of main rainy season and length of the growing period (LGP) i.e. the period from short rain season to the end of main season) and the NDVI. Different methodologies were adopted to achieve the objective of the study.

2.3.1. Correlation Analysis of Rainfall Data

First and foremost, the study evaluated the performance of the reanalysis data from ERA5 and CHIRPS with the ground-based observational data for two sites in NEB. The correlation was systematically calculated between two different datasets employing Pearson Correlation Coefficient (PCC) which takes values within the range [−1, 1] and is calculated with:

$r_{xy}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{\sqrt{{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left(x_i-\bar{x}\right)}^2{\left(y_i-\bar{y}\right)}^2}},$ (1)

where x, y, and n are independent variable (observed monthly rainfall data), dependent variable (predicted monthly rainfall, ERA and CHIRPS), and the number of observations times, respectively.

In addition, to understand the role of climatic factors on vegetation dynamics, a PCC analysis between MODIS NDVI and climatic variables (precipitation, temperature, onset of growing period, end of growing period, length of growing period and evapotranspiration) was computed using the commonly used approach, as shown in Equation (1), where $r_{xy}$ is the PCC of x and y variables, $y_i$ is the dependent variable representing mean annual or seasonal NDVI of the i^th year, and $x_i$ is the independent variable representing climate variable of the i^th year. $\bar{y}$ represents the average NDVI for the study period (18 years) and $\bar{x}$ represents the mean of climate factor of interest for the study period. The PCC is evaluated using F-test and a p-value less than 0.05 was considered statistically significant.

2.3.2. Trend Analysis

A number of tests are available to detect and estimate trends for temperature and rainfall data which often indicate the long-term change pattern. For instance, the nonparametric Mann-Kendall (MK) test (Mann, 1945; Kendall, 1975) has been proposed as an appropriate statistical tool to assess the data trends in time series of environmental variables. It is independent of the data distribution and consists of comparing each value of the time series with the other values remaining in the sequential order. The MK statistic S is defined as follows (Salmi et al., 2002):

$S={\displaystyle \underset{i=j+1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}}\text{sign}\left(x_i-x_j\right)\mathrm{,}$ (2)

where $x_i\mathrm{,}{x}_j$ are the values of the sequential data point ranked from $i=1,2,\cdots ,n-1;j=1,2,\cdots ,n$ ; n is is the length of the time series; $\text{sign}\left(x_i-x_j\right)$ is the function and is given as:

$\text{sign}\left(x_i-x_j\right)=\{\begin{array}{l}1\text{if}\left(x_i-x_j\right)>0\\ 0\text{if}\left(x_i-x_j\right)=0\\ -1\text{if}\left(x_i-x_j\right)<0\end{array}$ (3)

The mean $E\left[S\right]$ and the variance $V\left[S\right]$ under the null hypothesis (i.e. there is no trend) are calculated as (Kendall, 1975);

$\begin{array}{l}E\left[S\right]=0\\ V\left[S\right]=\frac{1}{8}\left\{n\left(n-1\right)\left(2n+5\right)-{\displaystyle \underset{i=1}{\overset{m}{\sum }}}{t}_i\left(t_i-1\right)\left(2t_i+5\right)\right\},\end{array}$ (4)

where $t_i$ is considered as the number of connections to the i^th value, and m the number of connected groups. The test statistics Z is calculated using $V\left[S\right]$ as follows:

$Z=(\begin{array}{ll}\frac{S-1}{\sqrt{V\left[S\right]}}\hfill & \text{if}S>0,\hfill \\ 0\hfill & \text{if}S=0,\hfill \\ \frac{S+1}{\sqrt{V\left[S\right]}}\hfill & \text{if}S<0,\hfill \end{array}$ (5)

The computed value of Z is then used to analyze the presence of a statistically significant trend. A positive Z value indicates a positive trend, while a negative Z value points to a negative trend.

The slope of precipitation trend was computed using Sen’s estimator (Sen, 1968). According to this method, the trend is assumed to be linear, i.e.

$f\left(t\right)=Qt+B,$ (6)

where f(t) is an increasing or decreasing function of time, with Q the trend’s slope and B the intercept (constant). The slope of each data pair $Q_i$ is calculated as:

$Q_i=\frac{x_j-x_k}{j-k},$ (7)

where $j>k$ and, if there is n number of $x_j$ in the time series, we get as many as $N=\frac{n\left(n-1\right)}{2}$ slope estimates of $Q_i$ . Then the values of $Q_i$ are ranked from small to large; the median of which is the Sen’s slope (Q):

$Q=(\begin{array}{ll}{Q}_{\left[\frac{n+1}{2}\right]}\hfill & \text{if}N\text{is}\text{odd}\hfill \\ \frac{1}{2}\left(Q_{\left[\frac{n}{2}\right]}+Q_{\left[\frac{n+2}{2}\right]}\right)\hfill & \text{if}N\text{is}\text{even}\hfill \end{array}$ (8)

These MK’s test and Sen’s Slope Estimator test, which have gained acceptance in hydro-climatology, was used for yearly precipitation data trend detection.

2.3.3. Variability Analysis

Temporal climate variability was determined using the coefficient of variation (CV) which is computed as:

$\text{CV}=\left(\frac{\text{SD}}{\text{Mean}}\right)\times \mathrm{100,}$ (9)

where CV represents the coefficient of variation; SD is the standard deviation of a station for the period of analysis and Mean is the mean rainfall for period of analysis at a given station. A CV is the ratio of the standard deviation to the mean of the rainfall at any given station (Egeru et al., 2014). Classifying the CV% has been the subject of interest in several areas, where researchers proposed different categories to rank it. Hare (1983) classified the CV as low for below 20% CV, high for CV in the range of 20% - 30% and more than 30% as very high variability. On the other hand, Vörösmarty et al. (2005) state that stations with CVs less than 30% show a moderate variability.

2.3.4. Precipitation Concentration Index

The Precipitation Concentration Index (PCI) as introduced by Oliver (1980) and further developed by de Luis et al. (2011) has been expressed as an indicator of rainfall concentration for annual scales. The PCI is used to examine 32-year rainfall data for the NEB area. The equations used to analyze the PCI are written as follows (de Luis et al., 2011; Oliver, 1980):

$\text{PCI}=\frac{{\displaystyle \underset{i=1}{\overset{12}{\sum }}}{P}_i^2}{{\left({\displaystyle \underset{i=1}{\overset{12}{\sum }}}{P}_i\right)}^2}\times 100,$ (10)

where $P_i$ is the rainfall amount of the i^th month. In addition, the PCI was calculated on a decadal scale for 1986-1995, 1996-2005 and 2006-2015, using mean monthly rainfall data splited into these decades respectively. Oliver (1980) and de Luis et al. (2011) classified PCI values below 10 as a uniform monthly rainfall distribution (low precipitation concentration); values between 11 and 20 as high concentrations of monthly rainfall distribution; and values of 21 and above as very high concentrations of monthly rainfall distribution (strong irregularity of precipitation distribution).

2.3.5. Growing Season Characteristics

The onset and cessation of rainfall strongly determine the agronomic processes in the study region. Yet, many definitions exist to define these key agroclimatic features.

Stern et al. (1982) defined the start of the rainy season as the first occurrence of at least a quantity of rainfall totaled over a number of consecutive days. This potential start has been shown false when a dry spell of “n” or more days in the next “m” days occurs. Accordingly, the earliest start of the growing season is the first occasion when the rainfall accumulated within a 3-day period exceeds 20 mm or more. Different authors have used a similar definition for assessing the start of the rainy season (Barron et al., 2003; Stern et al., 2006; Kassie et al., 2014; Nanja, 2010). For instance, Edoga (2007) defines the start of rain when the monthly rainfall is at least 60 mm. Dorward et al. (2015) define it as when 20 mm of rain precipitates within the next thirty days and as at least 20 mm is not followed by a period of more than 10 consecutive dry days.

For this study, based on the definition of the onset of the rains provided by Dorward et al. (2015), two scenarios were considered to calculate the onset of the rainy season: “start” without considering any “dry spells”, and “stdry” being “start with dry spells”. The “start” was defined as the first occasion with more than 20 mm in 2 consecutive days after 1st September while “stdry” was defined as the first occasion with more than 20 mm in 2 days after 1st September where no dry spells occur that exceed 9 days after the first rain. The growing period is referred to the period from onset of the short rainy season to the end of the main rainy season.

The end of the main rainy season was defined as the moment when starting from first May, the soil water is depleted (i.e. when the soil water balance reaches zero) by illustrating here a definition based on a simple water balance equation (Stern et al., 1982):

Today’s Water Balance

= Yesterday’s Water Balance + Precipitation − Evaporation; (11)

where evaporation is estimated to 5 mm. Thus, based on the annual values of the onset and cessation of the rainy season, we calculated the length of the growing period.

Climatic data analysis has traditionally relied on monthly rainfall to understand annual or seasonal rainfall patterns. Different methods are used to study a one-dimensional data set. For instance, statistical boxplots have been shown useful in establishing the seasonal and interannual variabilities of climatic events. Boxplots include information on estimated values, their location (mean or median), scale (interquartile range), and asymmetry (difference between quartile and median).

2.3.6. Dry Spell Analysis

A simple Markov chain model was used to model all the relevant information from the rainfall data regarding the probabilities of rain and to the rainfall amounts. A Markov chain consists of describing the chance of rain dependently to the previous day status i.e. the chance that a dry spell continues after a dry day and also the chance that the rain spell continues after a rainy day (Stern et al., 2006; Stern & Cooper, 2011).

Let ${\left\{X_n\right\}}_{n\ge 0}$ be a stochastic process that takes a finite number of possible values and S be a set of possible states. If $X_n=i$ with $i\in S$ , the process is said to be in state i at time n. We suppose that, whenever the process is in state i, there is a fixed probability, $P_{ij}$ , that it will be in state j (Ross , 2014 ). We have the transition matrix P defined as follows:

$P={\left(P_{ij}\right)}_{\left(i\mathrm{,}j\in S\right)}\mathrm{.}$ (12)

Consider 1 and 0 as the conditions of the day where a day is considered a rain (wet) day if the rainfall received is at least 0.85 mm, and a dry day otherwise. The random variables $X_n$ take the values 0 and 1 over T days in period of study. Thus, with $S=\left\{\mathrm{0,1}\right\}$ :

$X_n=(\begin{array}{ll}0\hfill & \text{if}\text{the}n\text{-th}\text{day}\text{is}\text{dry,}\hfill \\ 1\hfill & \text{if}\text{the}n\text{-th}\text{day}\text{is}\text{wet,}\text{with}n=0,1,2,\cdots ,T\mathrm{.}\hfill \end{array}$ (13)

Let us assume that, the probability that it will rain today depends on the state of the previous day only, that is:

$P\left(X_{n+1}=j|X_n=i,X_{n-1}=i_{n-1},\cdots ,X_0=i_0\right)=P\left(X_{n+1}=j|X_n=i\right)=P_{ij}\left(n\right),$ (14)

where $i_0\mathrm{,}{i}_1\mathrm{,}\cdots \mathrm{,}{i}_{n-1}\mathrm{,}j\in \left\{\mathrm{0,1}\right\}$ and $P_{ij}\left(n\right)$ is the transition probability with n steps (Ross, 2014). In other words, we assume that the probability of wetness on any day only depends on whether the previous day was wet or dry. This first order Markov chain stochastic process stipulates that, given the event on the previous day, the probability of wetness is assumed independent of further preceding days (Ross, 2014).

The probability of a rainy day was assessed by computing for each day of the year the frequency of rainfall greater than 0.85 mm. The frequency of dry spells was calculated for two periods of the rainy period: the short rainy season and the main rainy season. To define the length of the dry spell period and its occurrence, for each period, the frequency of 3, 5, 8, and 11 consecutive days with little (less than 0.85 mm/day) intensity or no rain was considered. It was calculated for each day of the rainy seasons by clustering the periods after the considered day, using the first order Markov chain model (Kassie et al., 2014).

2.3.7. Remote Sensing of Vegetation Cover

The study of vegetation using satellite imagery is made possible by the spectral characteristics of the plants. Indeed, plants have particular behaviors with regard to reflectance in the electromagnetic spectrum. In the visible range (400 - 700 nm), most of the radiation is absorbed by leaf pigments (chlorophyll, carotene, among others). The higher the level of photosynthesis, the lower the reflectance (Harhouz, 2012).

In the near-infrared range (700 - 1300 nm), the leaf pigments, as well as the cellulose, are transparent. The received radiation is, therefore, either reflected or transmitted. By increasing the wavelengths from the visible to the near-infrared, the reflectance is changed from very low to near 40%. In the mid-infrared range, the reflectance of the plants is mainly affected by their moisture content (Girard & Girard, 2010). The red (R) and near-infrared (NIR) bands show a low dependence on atmospheric conditions. The first coincides with strong absorption of radiation, the second offers high reflectance. Thus, it is possible to construct indices from a simple or complex combination of spectral bands. The combinations of these bands have good discriminatory capabilities and reveal particular properties of plants, for example, their chlorophyll content. The difference between R and NIR has led to the production of numerous vegetation indices highlighting the value of this difference in order to measure the photosynthetic activity of the plant (Girard & Girard, 2010). The mostly used indices are presented in Table 2.

These indices are highly correlated with the photosynthetic activity of biomass, chlorophyll abundance and uptake of energy by plants (Myneni et al., 1995). Despite the development of several vegetation indices, the NDVI remains the most widely used. This index has shown a strong correlation with vegetation parameters such as the vegetation vitality or plant biomass and is considered a good vegetation discriminator (Kennedy, 1989). Thus, NDVI is used to analyze seasonal changes in vegetation because its sensitivity is interesting especially for low-dense canopies (Caloz & Collet, 2001).

To assess the evolution of vegetation growth throughout the period of study (i.e. from February 2000 to December 2017), the MODIS 16-day composite NDVI values were averaged at civil year scale (NDVI year) and seasonal scale (SOND average: NDVI short, and MAM average: NDVI main). Key climatic variables linked to NDVI include start of short rainy season, end of the main rainy season, length of the growing period, mean annual air temperature, annual precipitation, and mean annual potential evapotranspiration. The relationship between the NDVI and climate variables was analyzed in two seasons (short rainy season and main rainy season) of plant growth separately. Thus, cumulative precipitation, average NDVI, average air temperature, and averaged evapotranspiration values were taken into account for these seasons.

Instat plus Version 3.37 software (Stern et al., 2006) was used to determine the onset, cessation and length of the growing season, the annual and seasonal rainfall totals, the annual seasonal number of rain days and the maximum of spell

Table 2. The most commonly used vegetation indices. R: Red channel; NIR: Near infrared channel; L: Soil adjustment factor; a: slope of the blank line in the R/NIR space; RB*: mean apparent reflectance of the red channel; NIR*: mean apparent reflectance of the near infrared channel (modified from Quense, 2011).

length while R (R Core Team, 2020) and python have been used for other analyses.

3. Results and Discussion

3.1. Rainfall Characteristics and Variability

Different methods are used to analyze the climatology of NEB. The annual cycle of rainfall and correlation between reanalysis and ground-based station rainfall data was analyzed to assess their performance. It was found that gridded data from ECMWF atmospheric reanalysis datasets and CHIRPS data are in agreement with the station data with a correlation of 0.77 and 0.84 in Kirundo and 0.81 and 0.88 in Muyinga for ERA5 and CHIRPS, respectively. Comparing these satellite products, CHIRPS performs better than ERA5. The CHIRPS dataset has been proven valid in East Africa (Dinku et al., 2018), and has been used in diverse works (Nkunzimana et al., 2019; Muthoni et al., 2019; Cattani et al., 2018; Belay et al., 2019).Figure 2 shows the correlation between the datasets used through two different stations which are Kirundo and Muyinga.

Figure 2. Correlation coefficient of monthly gridded of two rainfall zones over NEB based on CHIRPS and ERA5 datasets for a period of 32 years (1986-2017).

Burundi has various agro-ecological regions with varying patterns of precipitation. However, all regions have a bimodal rainfall pattern with peaks in April (Mid-February-May, main season) and November (September-December, short rain season). The bimodal pattern of rainfall for all data datasets is given in Figure 3.

The consecutive months JJA are characterized by low rainfall at both Kirundo and Muyinga stations and correspond to the dry season. A minimum of rainfall in January and February is observed, that characterizes a short dry season. The differences in mean annual precipitation in the studied area may be explained by their geographical positions with respect to the altitude. The highlands are ranked in the climatic zone with three dry months while the lowlands account four to six dry months (Lawin et al., 2018).

The long-term annual rainfall total varied from 496 mm to 1451.6 mm in Kirundo station. The main growing season received, on average, 51.3% of the annual rainfall total while the short season contributed to 48% of the total rain received by that region. For Muyinga, the range of 769.2 mm to 1580.9 mm is recorded in annual rainfall total with 49% for the main growing season and 39.1% for the short season, the rest falling in the dry season. The Muyinga station has recorded higher seasonal and annual rainfall than Kirundo. Rainfall total may not be sufficient to help agricultural planners for rain-fed crop production.

From Table 3, it can be noted that both stations fall under low variability (CV < 20%) in terms of their annual total rainfall and main season while they fall under high variability (CV > 20%) in terms of short rain season (Hare, 1983). The short rain season exhibits a high risk of rain-fed agriculture production failure. Similar results have been provided by Baramburiye et al. (2012) for NEB, who suggested to introduce drought-tolerant agricultural technologies to adapt to climate change.

The inter-annual variabilities are due to the annual movement of the Inter-Tropical Convergence Zone (ITCZ) which varies predictably throughout the year (Misiani et al., 2019). Variations in the ITCZ location significantly affect

Figure 3. Mean annual cycle of rainfall (mm) over the NEB for station data, ERA5 data and CHIRPS data for the period (1986-2017) and annual cycle of composite NDVI over the region for the period 2000-2017.

Table 3. Rainfall (mm) variability characteristics of two selected stations in the NEB, 1986-2017.

rainfall in many equatorial countries and lead to wet and dry seasons in the tropics. The ITCZ is less mobile over the ocean. Yet, the ITCZ migrates latitudinally on a seasonal basis (Misiani et al., 2019). In July, when the sun is over the Tropic of Cancer, the ITCZ reaches its northernmost position at about 15˚N latitude; in January it reaches 5˚S latitude when the sun is over the Tropic of Capricorn. The East African annual alteration of wet and dry seasons is mostly caused by the ITCZ shifting. When the ITCZ is to the north of the equator, the south westerly winds prevail over Burundi producing the dry-season conditions. Vegetation growth is generally restricted, grasses dry and leaves fall from trees due to reduced humidity. The ITCZ moves back into the Southern Hemisphere in September, bringing rain during the wet season.

The PCI values also indicated that in both stations the values were more than 10%, indicating the presence of seasonality in rainfall distribution. On an annual basis, the result of the PCI calculated for Kirundo showed that the PCI had the lowest value for 1989 (9.63%) and the highest for 2016 (16.21%) with an average of 13.3%. For Muyinga, an average of 13.34% is observed with a minimum of 10.81% in 2011 and a maximum of 17.16% in 1999. To investigate further the features of PCI, values were calculated based on decadal scale, where the decadal PCI for 1986-1995, 1996-2005, 2006-2015 falls into the category of moderate precipitation concentrations, giving an average value of 13.09, 12.9 and 13 respectively in Kirundo and 13.19, 13.93 and 12.9 respectively for Muyinga.

3.2. Annual and Seasonal Rainfall Trends in NEB

Changes and variabilities in annual and seasonal rainfall totals were assessed and are shown in Figure 4. The mean annual precipitation in the two locations was

Figure 4. Trends of annual and seasonal rainfall amount at Kirundo (left) and Muyinga (right) stations in NEB.

calculated, 1053.6 mm for Kirundo station and 1092.4 mm for Muyinga station. Using the Mann Kendal test and Sen’s slope for trend and its magnitude assessment, it is indicated that there is no-significant (p > 0.05) trends observed in annual rainfall, short season and main season in both stations, during the recorded period (Table 4). For the studied period (1986-2017), annual rainfall and main rainy season have decreased insignificantly (p > 0.05) over time at Kirundo and Muyinga with an average drop of 5.8 mm/year and 2.48 mm/year, respectively at Muyinga station and 0.6 mm/year and 0.71mm/year respectively at Kirundo. As for the short season, an average drop of 0.12 mm/year is observed in Muyinga while an average rise of 2.36 mm/year is revealed in Kirundo for the period (1986-2017). There is no evidence observed of trend and climate change (all changes are no-significant at 95%), contrariwise a significant variability.

Although there are no significant trends in rainfall patterns, the high inter-annual variability and seasonal variation by season represent a challenge for rain-fed agriculture. Different studies have suggested that the extent and temporal distribution of rainfall is the most important determinant of year-to-year fluctuations in agricultural production in Burundi (Baramburiye et al., 2012), and has had a considerable influence on the national economy and food production over the past three decades (Lawin et al., 2018; Thornton et al., 2019).

3.3. Characteristics of the Growing Period

3.3.1. Onset, Cessation and Length of Growing Period

The onset dates of the rainy season vary according to the regions and are conceived in different ways by farmers, and according to their soils. Studies carried have shown that, besides the factors that challenge the agriculture productivity in rural areas (soil fertility, availability of labor, seed, and fertilizer), the availability of information on climate also influence the decision making of smallholder farming (Nanja, 2010; Beyer et al., 2016; Cooper et al., 2008; Calzadilla et al., 2013). However, the accessibility to this information is limited.

Table 4. Annual and seasonal rainfall trend (mm/yr) of selected stations in the NEB.

Onset dates are given in Figure 5. They range from September, 2nd to November, 1st with mean on September, 23th for “start” while they range from September 2 to November 30 with mean October 8 for “stdry”, for Kirundo. The earliest onset day as well for “start” and “stdry” occurred in 2004 while the latest onset day occurred in 2017 and 2000. For Muyinga station, the earliest onset date was September 3 and the latest date was November 26, for as well “start” as “stdry”. The standard deviation of onset was 16.3, 19 days respectively for “start” and “stdry” in Kirundo and 21.5 and 20.3 days for “start” and “stdry” in Muyinga. High standard deviations of onset date of seasons lead to a misunderstanding of climate dynamics and, as a result, decisions about planting and related activities will be made at high risk (Ademe et al., 2019). The onset of the growing season may be associated with the northward/southward movement of the ITCZ and the Congo air mass which are subjected also to variability in year scale (Mugalavai et al., 2008; Ininda, 1994).

Figure 5 shows that the start and stdry patterns for Kirundo and Muyinga are different for some years. The start pattern is lower than that of stdry due to the reason that the dry days occur later after start conditions when 1st September is taken as the first day. However, the two definitions led to same output from 1986 to 2006. Results revealed that 17 and 13 over 32 years (about 53% and 40.6%) experienced above normal starts of rain (i.e. 8th of October) in Kirundo and Muyinga, respectively. These percentages show that the delay in the start of the rainy season can be quite significant and should not be ignored. Therefore, these years are years in which risks of planting failures were high and farmers would have had to replant. The two definitions coincide from 1986 to 2001 in Muyinga and vary differently for other years while they gave different values for only three years (2007, 2014 and 2016) in Kirundo. On running the MK test and Sen’s slope for trend and its magnitude assessment, the results in Table 5 were obtained. The onset dates of rainy period tend to increase slightly (positive Z), meaning that the delays in the latter continue to be observed over time for both stations, for the period of study. However, the cessation dates of the rainy period have decreased as well as the length of the growing period.

3.3.2. Growing Period Characteristics

The relationship between growing period characteristics is analyzed. The mean length of the growing period, as shown in Table 6, is the same for both stations while the mean onset date differs for one week. The probability of occurrence of early, normal, and late onset and cessation is also given in Table 6. The data show

Figure 5. Start date of rain based on two definitions: “start” and “stdry” (date of year after 1^st September) against time (year) for Kirundo (top) and Muyinga (bottom) stations.

Table 5. Onset, cessation and length of growing period trend and variability for two selected stations.

Table 6. Mean date of onset, cessation, duration, and frequency of occurrence of normal, early, and late onset short rainy season and cessation of the main rainy season in two locations of the Northeast Burundi.

that the probability of early onset and early cessation is similar for two stations(19% and 22%) while Kirundo has a lower probability of late onset and a higher probability of late cessation than Muyinga. The correlation analysis between the onset and cessation of the growing period indicated that there is no significant relationship at both stations. A stronger significant (p < 0.001) inverse relationship between onset and length/duration of the growing season is revealed at both Kirundo and Muyinga stations (Figure 6). The presence of a relationship between the onset or cessation dates, on the one hand, and between the length of the season is important for the planning of activities in the season. This relation is consistent with the analysis carried out by Sivakumar (1988) on long-term daily rainfall data for 58 locations in the southern Sahelian and Sudanian climatic zones of West Africa. Oladipo and Kyari (1993), when investigating the fluctuation in the onset, cessation, and length of the growing season in Northern Nigeria, also reported the sensitivity of the growing season to the onset of rainfall.

Figure 6. Relationship between the onset, cessation, length of the rainy season and annual precipitation in two locations in NEB.

3.3.3. Patterns of Dry Spells during the Growing Period

To estimate the risk of a long dry spell after sowing, the probability of dry spell lengths of 3 days (sp3), 5 days (sp5), 8 days (sp8) and 11 days (sp11) that occurred during the short rain and main rain growing seasons was estimated based on the first-order Markov Chain Model (Figure 7 and Figure 8). The analysis revealed that dry spell lengths were more prevalent in Kirundo than in Muyinga. The probability of short dry spells (3 - 5 days) is high and varies from 50% to more than 95% in the short rain season for both stations. On the other hand, the risk of long dry spells (8 - 11 days) is relatively low between day of year (DOY) 260 to 340 and 60 to 120 in short rain season and main rain season respectively. In the short rain season, the risk of short dry spells increases fast after DOY 330 in Kirundo, and DOY 320 in Muyinga while a fast increase of the risk of long dry spells is observed after DOY 340 in Muyinga and Kirundo.

As the months of January and February have high probability of getting dry spell of 3 days, they are considered as the small dry season (Nkunzimana et al., 2019) and Mid-February is taken as the average start date of the main rainy season. There is high probability (more than 90%) of 3 days dry spell, more than 50% of 5 dry spell in the second half of February and 100% 3 days dry spells after April 27. The results showed that if a crop could withstand a 10 day dry spell (but no longer) then waiting until about day 118 (27th April) before sowing, would give an estimated risk of just 20%, i.e. one year in five, that the crop would fail. If the crop could not withstand a dry spell of more than 7 days then the earliest sowing should be about 10 days later, i.e. not before 7th May.

Dry spells are important precipitation characteristics that drive the occurrence and duration of intra-seasonal droughts. Long dry spells cause significant yield loss if they coincide especially with drought-sensitive growth stages such as flowering and grain filling stages (Stern & Coe, 1984).

3.4. Changes in NDVI and Climatic Variables

The NDVI has been proven to be a suitable index to characterize the regional

Figure 7. Probability of dry spell occurrence of 3, 5 8 and 11 days of Muyinga station: short rain season (left) and main rain season (right).

Figure 8. Probability of dry spell occurrence of 3, 5 8 and 11 days of Kirundo station: short rain season (left) and main rain season (right).

distribution of vegetation and the seasonal dynamics of phenology (Kalisa et al., 2019; Ndayisaba et al., 2017). It is the most commonly used remote sensing index for assessing changes in vegetation conditions. Hence, changes in vegetation growth and vigour in NEB were assessed using NDVI. The spatial distribution of annual average NDVI in NEB from 2000 to 2017 is presented. The averaged NDVI’s trends, as well as its relationship with climate factors (precipitation, air temperature, and the potential evapotranspiration), are analyzed at annual and seasonal scales, to further investigate how and why vegetation occurred under water and energy fluctuations.

3.4.1. NDVI Dynamics over the Region

The trends over the 18 years when NDVI data were available were analyzed over the NEB.Figure 9 shows the interannual variation of annual mean NDVI from 2000 to 2017. Vegetation vigour along the period of study has generally been characterized by a slightly increasing trend (not significant). This reflects the improved agricultural and land-use practices promoted in recent decades. The minimum of 0.475 and a maximum of 0.591 are observed in 2000 and 2007 respectively. The linear regression model was used to assess the statistically insignificant and positive trend of 85 × 10−5 at p > 0.05. Yet, a breakpoint could be

Figure 9. Interannual variation of mean NDVI from 2000 to 2017.

detected using the piecewise linear regression model at the year 2011. Before the NDVI breakpoint, results show that the annual NDVI has increased at 40 × 10⁻⁴/yr and decreased at −12 × 10⁻³/yr after the breakpoint. In general, NDVI has been fluctuating. Figure 10 illustrates the spatial inter-annual variations of vegetation. The vegetation cover of the study area has been maintained at an average level except in 2000 where it is very low. However, areas of severe degradation extend along the northern lakes bordering Rwanda. Those undergoing improvement are mainly found in the northern depressions at Muyinga with altitudes varying between 1325 and 1430 m and in the Bugesera basin. This area is made up of lowlands interspersed by a group of hills rising to around 1550 m.

The southeasterly trade winds of the Indian Ocean are the main drivers of the climate of East Africa in general and the study area in particular and cause long (short) rains during MAM (SOND). The Intercontinental Convergence Zone (ICZ) has a biannual convergence with a seasonal inversion of the trade winds. The recession of the Indian monsoon system over East Africa limits the convection and southward flow arising from the northeastern highlands, resulting in southwestward moisture flow (Kalisa et al., 2019). As stated by Kalisa et al. (2019), most pluviometric irregularities involving flood and drought occurrences over the region are effects of El Niño Southern oscillation (ENSO) impacts. Yet, documented reports have shown that ENSO peak occurs between October and March, a period coinciding with rainy seasons in East Africa. It is worth recalling that the hot events of ENSO are contributing to increasing rainfall while the cold events lead to droughts affecting vegetation growth. To quantify the deficit (droughts) or excess (floods) of rainfall at different scales, the SPI was developed and shown in Figure 11. Thus, from Figure 10, it’s obvious that the period 2000 had a deficiency in Vegetation (low NDVI values) and coincides with a very strong 1998-2001 La Niña episode and the period 2007 with high NDVI coincides 2006-2007 moderate El Niño event (Null, 2017).

Figure 10. Spatiotemporal variability of NDVI year at the NEB for the period from 2000 to 2017. The regions with greener colors indicate higher NDVI values, whereas the brownish colors indicate low NDVI values.

Figure 11. Standarised precipitation index for twelve months time scales for the period from 1986 to 2017 in NEB based on CHIRPS rainfall dataset.

Although the variation trends of NDVI and climatic parameters were not significant (NDVI increase, water factors decrease and air temperature parameter increase), different studies stated that precipitation was the main factor determining vegetation growth (Gebru et al., 2020; Ndayisaba et al., 2017). During years with high precipitation, the vegetation becomes green while in years with low precipitation, the growth slows down and the vegetation becomes sparse. Excessive air temperature could limit vegetation growth because the increased air temperature leads to increased transpiration of vegetation, resulting in a lack of water available for vegetation growth, particularly in water-stressed areas (Li et al., 2003). In the study area, the mean NDVI first continued to increase until 2011/2012, before showing a downward trend until 2017. This could be due to drying out, through decreased precipitation and increased temperature. Although there was a decline in the mid-2011s after the breakpoint, the overall trend in NDVI remains positive as precipitation has continued to decline.

3.4.2. Relationship between NDVI and Climate Parameters at a Yearly Scale

NDVI and climate parameters of water and energy were normalized using the normalization method and their relationships were analyzed at a yearly scale. Changes in climate factors have influences on NDVI, and the change patterns of NDVI are for instance similar to those of precipitation. There is an increase in NDVI when precipitation is increased and decreased with decreased precipitation. Subsequently, the correlations between NDVI and meteorological parameters were analyzed quantitatively using the Pearson correlation coefficient and the partial correlation coefficient. The results showed that only precipitation had positive (not significant) correlations with NDVI (p > 0.05), and the other two parameters had negative correlations with NDVI (not significant) (Table 7).

Table 7. Pearson correlation coefficients and partial correlation coefficients between NDVI and meteorological parameters at a yearly scale in NEB.

Note: * represents $p<0.05$ significance; $R_{\text{NDVI-T}}$ , $R_{\text{NDVI-P}}$ and $R_{\text{NDVI-E}}$ represent correlation coefficients between NDVI and air temperature (˚C), NDVI and precipitation (mm), and NDVI and potentail evapotranspiration (kg/m²), respectively; and $P{R}_{\text{NDVI-T}}$ , $P{R}_{\text{NDVI-P}}$ , and $P{R}_{\text{NDVI-E}}$ represent partial correlation coefficients between NDVI and air temperature (˚C), NDVI and precipitation (mm), and NDVI and potentail evapotranspiration (kg/m²), respectively.

Thereupon, the partial correlation between air temperature and NDVI is revealed significant. These results do not show an evident influence of annual cumulative precipitation on vegetation growth in the region, nevertheless, seasonal fluctuations may be the main factor of variation in NDVI especially during the delays in the onset of rain observed in September.

3.4.3. Inter-Annual Changes in Seasonal Climate Dynamics and Their Relationship with Vegetation NDVI

The averaged NDVI for the period corresponding to the short rainy season and the main rainy season has been found to be year-to-year variable with no obvious trend over time. The coefficients of variation were respectively 6% and 3% for MAN and SOND.

Table 8 shows that the partial correlation between NDVI and meteorological parameters from high to low was as follows: air temperature > evapotranspiration > precipitation in the short rainy season; and evapotranspiration > precipitation > air temperature in main rain season. The difference in water and energy determinants at large time scales (e.g. annual scale) and small time scales (e.g. seasonal) suggests that the correlation between NDVI and meteorological parameters is scale-dependent.

Figure 12 illustrates the response of NDVI to climate parameters at yearly scale and in two seasons of plant growth separately. A significant negative correlation is revealed between NDVI year and the start of short rainy season. Also, NDVI short showed a significant negative correlation with SOND average air temperature and evapotranspiration. There was a tendency towards a positive correlation (significant at 90%; p < 0.10) between NDVI short and SOND cumulative precipitation. On the contrary, these variables have shown no significant correlation during MAN growth period. It has therefore been shown that the start date of the rainy season is the key variable to which the vegetation growth is mainly related. At the beginning of the growing period, the delays in the onset of rainfall, in combination with the increase in temperature and evapotranspiration, also plays a remarkable role in vegetation dynamics. Many farmers are very conscious of the fact that early onset rainy seasons are generally better for crop production than late onset rainy seasons (Sivakumar, 1988). The

Table 8. Pearson correlation coefficients and partial correlation coefficients between NDVI and meteorological parameters at a monthly scale in two growing seasons in NEB.

Note: * represents p-value < 0.05 and ** represents p-value < 0.01 significance; $R_{\text{NDVI-T}}$ , $R_{\text{NDVI-P}}$ and $R_{\text{NDVI-E}}$ represent correlation coefficients between NDVI and air temperature (˚C), NDVI and precipitation (mm), and NDVI and potential evapotranspiration (kg/m2), respectively; and $P{R}_{\text{NDVI-T}}$ , $P{R}_{\text{NDVI-P}}$ , and $P{R}_{\text{NDVI-E}}$ represent partial correlation coefficients between NDVI and air temperature (˚C), NDVI and precipitation (mm), and NDVI and potential evapotranspiration (kg/m2), respectively.

uncertainty about the relationship between the onset and the end of the rainy season is a challenge to the popular belief that the late onset of the rainy season is compensated for by the late end of the rainy season or that the rainy seasons become shorter due to a late onset and early end (Traore et al., 2013). The north-eastern region of Burundi, mainly consisting of lowlands areas, is characterized by a delay in the starting of rains compared to other regions and the dry season lasts for 4 to 6 months (Lawin et al., 2018).

As indicated in Figure 3, on average, the growing season starts well with the second half of September and ends by the second week of May with a breakdown in February. The accumulated precipitation maintains the vegetation alive until the second week of June and progressively decreases until the complete dryness in August. With the possible onset of rains in mid-September, the cover would not recover well until the second week of October. So the delayed onset of the rains will impact vegetation production during the first cropping season in north-eastern Burundi. This has already been evident in the communes neighboring the northern lakes. The description of NDVI’s seasonal changes reflects the mean dynamics of NDVI during the growing season. However, from an interannual perspective, the NDVI values, i.e. vegetation productivity, vary from one year to another. It is notable that there is a high inter-annual variability with unproductive seasons as in 2000-2001 and more favorable years as in 2006-2007.

4. Conclusion

This study analyzed the temporal variations of vegetation and their relation with meteorological parameters in NEB, a region presenting considerable climatic irregularities that have caused food shortages over recent decades. The 17-year remote sensing data of NDVI from the MODIS sensors over the surface vegetation of NEB were assessed in detail to study significant patterns in NDVI changes. The results showed that the spatiotemporal vegetation productivity and variation over the region of study vary across areas and time periods. The linear regression model and the piecewise regression model are robust to detect stepwise and breakpoints in mean NDVI and have been used to assess the statistical

Figure 12. Seasonal response of NDVI to water and temperature parameters in NEB. (top) represent relationships between NDVI year and onset dates, end of growing season and length of the growing season respectively; (middle) represent relationships between NDVI short and air temperature, precipitation and evapotranspiration respectively, during the short rainy season; and (bottom) represent relationships between NDVI main and air temperature, precipitation, and evapotranspiration respectively, in the main rainy season.

significance of trends. Results show an insignificant positive trend of 85 × 10⁻⁵ at p > 0.05. The breakpoint detected using the piecewise linear regression model corresponds to the year 2011. Before the NDVI breakpoint, results show that the annual NDVI has increased at 40 × 10⁻⁴/yr and decreased at −12 × 10⁻³/yr after the breakpoint. NDVI showed a positive (negative) relationship with precipitation (temperature and evapotranspiration) during the short rainy season while it shows reverse for the long rainy season.

Regarding the phenology metrics, results showed that, on average, the growing season starts well with the second half of September and ends by the second week of May with a breakdown in February. The accumulated precipitation maintains the vegetation alive until the second week of June and progressively decreases until the complete dryness in August. With the possible onset of rains in mid-September, the vegetation cover would not recover well until the second week of October. NDVI showed lags with precipitation on both long and short rainy seasons, which explains the fact that if the rains stop, vegetation activity persists before returning to its dormant state. It is revealed that the delayed onset of the rains has a significant impact on vegetation production during the first cropping season in north-eastern Burundi. Also, some relationship between NDVI and ENSO during most El Niño and La Niña years are observed, particularly in 1998/2001 and 2006/2007. The period 2000 had a deficiency in Vegetation (low NDVI values) and coincides with a very strong 1998-2001 La Niña episode and the period 2007 with high NDVI coincides 2006-2007 moderate El Niño event.

While other studies on vegetation response to climatic conditions stipulate a strong relationship between vegetation and precipitation, the regression analysis carried out has revealed that, in the study area, the major problems in vegetation production are mainly due to delayed rains in September and dry spell days occurring during the growing period. However, although this study contributed to the understanding of the complex responses of vegetation to climate conditions, it could not distinguish the rate of contribution of each meteorological factor. Other statistical methods, such as Random Forest (RF), wavelet analysis, and Causal effect Network (CEN) can be used in future studies.

Acknowledgements

This work was produced as part of a Doctoral Research of P. Batungwanayo, supported by the Doctoral School of University of Burundi and the Belgian cooperation through ARES (Académie de recherche et d’enseignement supérieur). M. Vanclooster contribution has been supported by the DGD-ARES-PSR program KLIMSEC. The authors are also gratefull to the Burundian National Meteorological Agency (IGEBU), the fifth generation of ECMWF atmospheric reanalyses and the climate hazard group for kindly providing the climate data. All other anonymous contributors whose comments helped to improve this paper are acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Abdulraheem, M., Charles, E., Owoade, O., & Onifade, A. (2017). An Investigation into the Effect of Rainfall Variability on Farming Activities in Rural Area. Asian Journal of Advances in Agricultural Research, 4, 1-7. https://doi.org/10.9734/AJAAR/2017/37872
[2]	Ademe, F., Kibret, K., Beyene, S., Getinet, M., & Mitikee, G. (2019). Rainfall Analysis for Rain-Fed Farming in the Great Rift Valley Basins of Ethiopia. Journal of Water and Climate Change, 11, 812-828. https://doi.org/10.2166/wcc.2019.242
[3]	Baramburiye, J., Kyotalimye, M., Thomas, T. S., & Waithaka, M. (2012). East African Agriculture and Climate Change: Burundi. In M. Waithaka, G. C. Nelson, T. S. Thomas, & M. Kyotalimye (Eds.), East African Agriculture and Climate Change (pp. 55-88). Washington DC: IFPRI.
[4]	Baret, F., Guyot, G., & Major, D. (1989). TSAVI: A Vegetation Index Which Minimizes Soil Brightness Effects on LAI and APAR Estimation. In 12th Canadian Symposium on Remote Sensing Geoscience and Remote Sensing Symposium (pp. 1355-1358). https://doi.org/10.1109/IGARSS.1989.576128
[5]	Barron, J., Rockström, J., Gichuki, F., & Hatibu, N. (2003). Dry Spell Analysis and Maize Yields for Two Semi-Arid Locations in East Africa. Agricultural and Forest Meteorology, 117, 23-37. https://doi.org/10.1016/S0168-1923(03)00037-6
[6]	Belay, A. S., Fenta, A. A., Yenehun, A., Nigate, F., Tilahun, S. A., Moges, M. M., Dessie, M., Adgo, E., Nyssen, J., Chen, M. et al. (2019). Evaluation and Application of Multi-Source Satellite Rainfall Product CHIRPS to Assess Spatio-Temporal Rainfall Variability on Data-Sparse Western Margins of Ethiopian Highlands. Remote Sensing, 11, 2688. https://doi.org/10.3390/rs11222688
[7]	Beyer, M., Wallner, M., Bahlmann, L., Thiemig, V., Dietrich, J., & Billib, M. (2016). Rainfall Characteristics and Their Implications for Rain-Fed Agriculture: A Case Study in the Upper Zambezi River Basin. Hydrological Sciences Journal, 61, 321-343. https://doi.org/10.1080/02626667.2014.983519
[8]	Bidou, J.-E., Ndayirukiye, S., Ndayishimiye, J.-P., & Sirven, P. (1991). Géographie du burundi.
[9]	Caloz, R., & Collet, C. (2001). Précis de télédétection—Volume 3: Traitements numériques d’images de télédétection. Québec: Presses de l’Université du Québec. https://doi.org/10.2307/j.ctv5j018b
[10]	Calzadilla, A., Zhu, T., Rehdanz, K., Tol, R. S., & Ringler, C. (2013). Economywide Impacts of Climate Change on Agriculture in Sub-Saharan Africa. Ecological Economics, 93, 150-165. https://doi.org/10.1016/j.ecolecon.2013.05.006
[11]	Cattani, E., Merino, A., Guijarro, J. A., & Levizzani, V. (2018). East Africa Rainfall Trends and Variability 1983-2015 Using Three Long-Term Satellite Products. Remote Sensing, 10, 931. https://doi.org/10.3390/rs10060931
[12]	Cooper, P., Dimes, J., Rao, K., Shapiro, B., Shiferaw, B., & Twomlow, S. (2008). Coping Better with Current Climatic Variability in the Rain-Fed Farming Systems of Sub-Saharan Africa: An Essential First Step in Adapting to Future Climate Change? Agriculture, Ecosystems & Environment, 126, 24-35. https://doi.org/10.1016/j.agee.2008.01.007
[13]	de Luis, M., Gonzalez-Hidalgo, J., Brunetti, M., & Longares, L. (2011). Precipitation Concentration Changes in Spain 1946-2005. Natural Hazards and Earth System Sciences, 11, 1259-1265. https://doi.org/10.5194/nhess-11-1259-2011
[14]	Dinku, T., Ceccato, P., Grover-Kopec, E., Lemma, M., Connor, S., & Ropelewski, C. (2007). Validation of Satellite Rainfall Products over East Africa’s Complex Topography. International Journal of Remote Sensing, 28, 1503-1526. https://doi.org/10.1080/01431160600954688
[15]	Dinku, T., Funk, C., Peterson, P., Maidment, R., Tadesse, T., Gadain, H., & Ceccato, P. (2018). Validation of the CHIRPS Satellite Rainfall Estimates over Eastern Africa. Quarterly Journal of the Royal Meteorological Society, 144, 292-312. https://doi.org/10.1002/qj.3244
[16]	Dorward, P., Clarkson, G., & Stern, R. (2015). Participatory Integrated Climate Services for Agriculture (PICSA): Field Manual. Reading: Walker Institute, University of Reading.
[17]	Edoga, R. N. (2007). Determination of Length of Growing Season in Samaru Using Different Potential Evapotranspiration Models. AU Journal of Technology, 11, 28-35.
[18]	Egeru, A., Osaliya, R., Macopiyo, L., Mburu, J., Wasonga, O., Barasa, B., Said, M., Aleper, D., & Majaliwa Mwanjalolo, G.-J. (2014). Assessing the Spatio-Temporal Climate Variability in Semi-Arid Karamoja Sub-Region in North-Eastern Uganda. International Journal of Environmental Studies, 71, 490-509. https://doi.org/10.1080/00207233.2014.919729
[19]	Field, C. B., Barros, V. R., Dokken, D. J., Mach, K. J., Mastrandrea, M. D., Bilir, T. E., Chatterjee, M., Ebi, K. L., Estrada, Y. O., Genova, R. C. et al. (2017). Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Cambridge: Cambridge University Press.
[20]	Funk, C., Peterson, P., Landsfeld, M., Pedreros, D., Verdin, J., Shukla, S., Husak, G., Rowland, J., Harrison, L., Hoell, A. et al. (2015). The Climate Hazards Infrared Precipitation with Stations—A New Environmental Record for Monitoring Extremes. Scientific Data, 2, Article No. 150066. https://doi.org/10.1038/sdata.2015.66
[21]	Gebru, B. M., Lee, W.-K., Khamzina, A., Wang, S. W., Cha, S., Song, C., & Lamchin, M. (2020). Spatiotemporal Multi-Index Analysis of Desertification in Dry Afromontane Forests of Northern Ethiopia. Environment, Development and Sustainability, 1-28. https://doi.org/10.1007/s10668-020-00587-3
[22]	Girard, M.-C., & Girard, C.-M. (2010). Traitement des données de télédétection—2e édition. Environnement et ressources naturelles. Dunod.
[23]	Hare, F. (1983). Climate and Desertification. A Revised Analysis. WCP-44. World Climate Program. Geneva: World Meteorological Organization/United Nations Environment Programme.
[24]	Harhouz, A. (2012). évaluation Multidate de la Qualité de l’Imageur d’Alsat-1. PhD Thesis, Bir El Djir: USTO.
[25]	Huete, A. R. (1988). A Soil-Adjusted Vegetation Index (SAVI). Remote Sensing of Environment, 25, 295-309. https://doi.org/10.1016/0034-4257(88)90106-X
[26]	Ininda, J. M. (1994). Numerical Simulation of the Influence of the Sea Surface Temperature Anomalies on the East African Seasonal: Rainfall. PhD Thesis, Nairobi: University of Nairobi.
[27]	Jordan, C. F. (1969). Derivation of Leaf-Area Index from Quality of Light on the Forest Floor. Ecology, 50, 663-666. https://doi.org/10.2307/1936256
[28]	Kalisa, W., Igbawua, T., Henchiri, M., Ali, S., Zhang, S., Bai, Y., & Zhang, J. (2019). Assessment of Climate Impact on Vegetation Dynamics over East Africa from 1982 to 2015. Scientific Reports, 9, Article No. 16865. https://doi.org/10.1038/s41598-019-53150-0
[29]	Kassie, B., Rötter, R., Hengsdijk, H., Asseng, S., Van Ittersum, M., Kahiluoto, H., & Van Keulen, H. (2014). Climate Variability and Change in the Central Rift Valley of Ethiopia: Challenges for Rainfed Crop Production. The Journal of Agricultural Science, 152, 58-74. https://doi.org/10.1017/S0021859612000986
[30]	Kaufman, Y. J., & Tanre, D. (1992). Atmospherically Resistant Vegetation Index (ARVI) for EOS-MODIS. IEEE Transactions on Geoscience and Remote Sensing, 30, 261-270. https://doi.org/10.1109/36.134076
[31]	Kendall, M. (1975). Rank Correlation Methods (4th ed.). San Francisco, CA: Charles Griffin.
[32]	Kennedy, P. (1989). Monitoring the Vegetation of Tunisian Grazing Lands Using the Normalized Difference Vegetation Index [AVHRR NDVI]. Ambio (Sweden), 18, 119-123.
[33]	Lawin, A. E., Manirakiza, C., & Lamboni, B. (2018). Trends and Changes Detection in Rainfall, Temperature and Wind Speed in Burundi. Journal of Water and Climate Change, 10, 852-870. https://doi.org/10.2166/wcc.2018.155
[34]	Li, Z., Yan, F., & Fan, X. (2003). The Variability of NDVI over Northwest China and Its Relation to Temperature and Precipitation. 2003 IEEE International Geoscience and Remote Sensing Symposium. Proceedings, Toulouse, 21-25 July 2003, 2275-2277.
[35]	Liersch, S., Rivas, R., & Fritzsche, K. (2014). Rapport sur le changement climatique au Burundi.
[36]	Mann, H. B. (1945). Nonparametric Tests against Trend. Econometrica, 13, 245-259. https://doi.org/10.2307/1907187
[37]	Manyeruke, C., Hamauswa, S., & Mhandara, L. (2013). The Effects of Climate Change and Variability on Food Security in Zimbabwe: A Socioeconomic and Political Analysis. International Journal of Humanities and Social Science, 3, 270-286.
[38]	McNally, A., Arsenault, K., Kumar, S., Shukla, S., Peterson, P., Wang, S., Funk, C., Peters-Lidard, C. D., & Verdin, J. P. (2017). A Land Data Assimilation System for Sub-Saharan Africa Food and Water Security Applications. Scientific Data, 4, Article No. 170012. https://doi.org/10.1038/sdata.2017.12
[39]	Measho, S., Chen, B., Trisurat, Y., Pellikka, P., Guo, L., Arunyawat, S., Tuankrua, V., Ogbazghi, W., & Yemane, T. (2019). Spatio-Temporal Analysis of Vegetation Dynamics as a Response to Climate Variability and Drought Patterns in the Semiarid Region, Eritrea. Remote Sensing, 11, 724. https://doi.org/10.3390/rs11060724
[40]	Ministry for Land Management and Environment (2007). National Adaptation Plan of Action (NAPA).
[41]	Misiani, Z., Lun, Y., Niu, S., Lü, J., & Zhang, L. (2019). A GIS Based Approach to Back Trajectory Analysis and Mass Concentration & Dispositions of Aerosols in Nairobi, Kenya. Journal of Geoscience and Environment Protection, 7, 122-139. https://doi.org/10.4236/gep.2019.72009
[42]	Mugalavai, E. M., Kipkorir, E. C., Raes, D., & Rao, M. S. (2008). Analysis of Rainfall Onset, Cessation and Length of Growing Season for Western Kenya. Agricultural and Forest Meteorology, 148, 1123-1135. https://doi.org/10.1016/j.agrformet.2008.02.013
[43]	Muthoni, F. K., Odongo, V. O., Ochieng, J., Mugalavai, E. M., Mourice, S. K., Hoesche- Zeledon, I., Mwila, M., & Bekunda, M. (2019). Long-Term Spatial-Temporal Trends and Variability of Rainfall over Eastern and Southern Africa. Theoretical and Applied Climatology, 137, 1869-1882. https://doi.org/10.1007/s00704-018-2712-1
[44]	Myneni, R. B., Hall, F. G., Sellers, P. J., & Marshak, A. L. (1995). The Interpretation of Spectral Vegetation Indexes. IEEE Transactions on Geoscience and Remote Sensing, 33, 481-486. https://doi.org/10.1109/TGRS.1995.8746029
[45]	Nanja, D. H. (2010). Dissemination of Climate Information to Small-Holder Farmers: A Case Study for Mujika Area, Zambia. Ph.D. Thesis, Bloemfontein: University of the Free State.
[46]	Ndayiragije, A., Mkezabahizi, D., Ndimubandi, J., & Kabogoye, F. (2017). A Scoping Study on Burundi’s Agricultural Production in a Changing Climate and the Supporting Policies. Kenya Institute for Public Policy Research and Analysis.
[47]	Ndayisaba, F., Guo, H., Isabwe, A., Bao, A., Nahayo, L., Khan, G., Kayiranga, A., Karamage, F., Muhire, E. N. et al. (2017). Inter-Annual Vegetation Changes in Response to Climate Variability in Rwanda. Journal of Environmental Protection, 8, 464-481. https://doi.org/10.4236/jep.2017.84033
[48]	Nicholson, S. E., Davenport, M. L., & Malo, A. R. (1990). A Comparison of the Vegetation Response to Rainfall in the Sahel and East Africa, Using Normalized Difference Vegetation Index from NOAA AVHRR. Climatic Change, 17, 209-241. https://doi.org/10.1007/BF00138369
[49]	Niragira, S., DHaese, M., DHaese, L., Ndimubandi, J., Desiere, S., & Buysse, J. (2015). Food for Survival: Diagnosing Crop Patterns to Secure Lower Threshold Food Security Levels in Farm Households of Burundi. Food and Nutrition Bulletin, 36, 196-210. https://doi.org/10.1177/0379572115587491
[50]	Nkunzimana, A., Bi, S., Jiang, T., Wu, W., & Abro, M. I. (2019). Spatiotemporal Variation of Rainfall and Occurrence of Extreme Events over Burundi during 1960 to 2010. Arabian Journal of Geosciences, 12, Article No. 176. https://doi.org/10.1007/s12517-019-4335-y
[51]	Null, J. (2017). El Niño and La Niña Years and Intensities Based on Oceanic Nino Index (ONI). https://ggweather.com/enso/oni.htm
[52]	Oladipo, E., & Kyari, J. (1993). Fluctuations in the Onset, Termination and Length of the Growing Season in Northern Nigeria. Theoretical and Applied Climatology, 47, 241-250. https://doi.org/10.1007/BF00866245
[53]	Oliver, J. E. (1980). Monthly Precipitation Distribution: A Comparative Index. The Professional Geographer, 32, 300-309. https://doi.org/10.1111/j.0033-0124.1980.00300.x
[54]	Pelkey, N., Stoner, C., & Caro, T. (2000). Vegetation in Tanzania: Assessing Long Term Trends and Effects of Protection Using Satellite Imagery. Biological Conservation, 94, 297-309. https://doi.org/10.1016/S0006-3207(99)00195-0
[55]	Pinty, B., & Verstraete, M. (1992). GEMI: A Non-Linear Index to Monitor Global Vegetation from Satellites. Vegetation, 101, 15-20. https://doi.org/10.1007/BF00031911
[56]	Quense, J. (2011). Changement climatique et dynamique de la végétation dans les Andes du Chili central, depuis le milieu du XXème siècle: L’exemple de la vallée de Yerba Loca. PhD Thesis, Saint-Martin-d’Hères: Université de Grenoble.
[57]	R Core Team (2020). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.
[58]	Richardson, A. J., & Wiegand, C. (1977). Distinguishing Vegetation from Soil Background Information. Photogrammetric Engineering and Remote Sensing, 43, 1541-1552.
[59]	Ross, S. M. (2014). Introduction to Probability Models. In Introduction to Probability Models (pp. 1-19). Cambridge, MA: Academic Press. https://doi.org/10.1016/B978-0-12-407948-9.00001-3
[60]	Rouse, J., Haas, R., Schell, J., & Deering, D. (1973). Monitoring Vegetation Systems in the Great Plains with ERTS. In Proceedings of the Earth Resources Technology Satellite Symposium (Vol. 1, pp. 309-317). Washington DC: NASA.
[61]	Salmi, T., Määttä, A., Anttila, P., Ruoho-Airola, T., & Amnell, T. (2002). Detecting Trends of Annual Values of Atmospheric Pollutants by the Mann-Kendall Test and Sen’s Solpe Estimates—The Excel Template Application MAKESENS. Helsinki: Finish Meteorological Institute.
[62]	Schreck III, C. J., & Semazzi, F. H. (2004). Variability of the Recent Climate of Eastern Africa. International Journal of Climatology:, 24, 681-701. https://doi.org/10.1002/joc.1019
[63]	Sen, P. K. (1968). Estimates of the Regression Coefficient Based on Kendall’s Tau. Journal of the American Statistical Association, 63, 1379-1389. https://doi.org/10.1080/01621459.1968.10480934
[64]	Sivakumar, M. (1988). Predicting Rainy Season Potential from the Onset of Rains in Southern Sahelian and Sudanian Climatic Zones of West Africa. Agricultural and Forest Meteorology, 42, 295-305. https://doi.org/10.1016/0168-1923(88)90039-1
[65]	Stern, R., & Coe, R. (1984). A Model Fitting Analysis of Daily Rainfall Data. Journal of the Royal Statistical Society. Series A (General), 147, 1-34. https://doi.org/10.2307/2981736
[66]	Stern, R., & Cooper, P. (2011). Assessing Climate Risk and Climate Change Using Rainfall Data—A Case Study from Zambia. Experimental Agriculture, 47, 241-266. https://doi.org/10.1017/S0014479711000081
[67]	Stern, R., Dennett, M., & Dale, I. (1982). Analysing Daily Rainfall Measurements to Give Agronomically Useful Results. I. Direct Methods. Experimental Agriculture, 18, 223-236. https://doi.org/10.1017/S001447970001379X
[68]	Stern, R., Rijks, D., Dale, I., & Knock, J. (2006). Instat Climatic Guide. Reading, UK: University of Reading.
[69]	Thornton, P., Loboguerrero, A., Campbell, B., Kavikumar, K., Mercado, L., & Shackleton, S. (2019). Rural Livelihoods, Food Security and Rural Transformation under Climate Change. Rotterdam and Washington DC: Global Commission on Adaptation.
[70]	Tirivangasi, H. M. (2018). Regional Disaster Risk Management Strategies for Food Security: Probing Southern African Development Community Channels for Influencing National Policy. Jàmbá: Journal of Disaster Risk Studies, 10, a468. https://doi.org/10.4102/jamba.v10i1.468
[71]	Traore, B., Corbeels, M., Van Wijk, M. T., Rufino, M. C., & Giller, K. E. (2013). Effects of Climate Variability and Climate Change on Crop Production in Southern Mali. European Journal of Agronomy, 49, 115-125. https://doi.org/10.1016/j.eja.2013.04.004
[72]	UNFCCC (2007). Climate Change: Impacts, Vulnerabilities and Adaptation in Developing Countries. Bonn: Climate Change Secretariat.
[73]	Vörösmarty, C. J., Douglas, E. M., Green, P. A., & Revenga, C. (2005). Geospatial Indicators of Emerging Water Stress: An Application to Africa. AMBIO, 34, 230-236. https://doi.org/10.1579/0044-7447-34.3.230
[74]	Yu, F., Price, K. P., Ellis, J., & Shi, P. (2003). Response of Seasonal Vegetation Development to Climatic Variations in Eastern Central Asia. Remote Sensing of Environment, 87, 42-54. https://doi.org/10.1016/S0034-4257(03)00144-5