Response of Seasonal Vegetation Dynamics to Climatic Constraints in Northeastern Burundi

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 spectrora-diometer (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, evapotranspiration, season characteristics length). show that the rate in vegetation productivity is persistently gradual 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.


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 Burun-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 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).

Area of Study
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   (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).

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. 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.

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.

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 da- 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 xy r is the PCC of x and y variables, i y is the dependent variable representing mean annual or seasonal NDVI of the i th year, and i x is the independent variable representing climate variable of the i th year. y represents the average NDVI for the study period (18 years) and 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.

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): where , i j x x are the values of the sequential data point ranked from 1, 2, , 1; 1, 2, , i n j n = − =   ; n is is the length of the time series; ( ) is the function and is given as: 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.
( ) 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 i Q is calculated as: where j k > and, if there is n number of j x in the time series, we get as many as ( ) Then the values of i Q are ranked from small to large; the median of which is the Sen's slope (Q): 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.

Variability Analysis
Temporal climate variability was determined using the coefficient of variation (CV) which is computed as: 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.

Precipitation Concentration Index
The Precipitation Concentration Index (PCI) as introduced by Oliver (1980) 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).

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; 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 Boxplots include information on estimated values, their location (mean or median), scale (interquartile range), and asymmetry (difference between quartile and median).

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).
n n X ≥ be a stochastic process that takes a finite number of possible values and S be a set of possible states. If n X i = with i 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, ij P , that it will be in state j (Ross, 2014). We have the transition matrix P defined as follows: 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 n X take the values 0 and 1 over T days in period of study.
Thus, with 0 if the -th day is dry, 1 if the -th day is wet, with 0,1, 2, , .
n n X n n T Let us assume that, the probability that it will rain today depends on the state of the previous day only, that is: ( ) ij P n 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). Journal of Geoscience and Environment Protection

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

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 agree-  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  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

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 P. Batungwanayo et al. 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)

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.  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 1 st 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. 8 th 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.

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    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.

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).

Changes in NDVI and Climatic Variables
The NDVI has been proven to be a suitable index to characterize the regional

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  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.

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). 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.

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 Journal of Geoscience and Environment Protection 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.

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 −5 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 −4 /yr and decreased at −12 × 10 −3 /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 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.