Spatio-Temporal Analysis of Drought Variability in Myanmar Based on the Standardized Precipitation Evapotranspiration Index (SPEI) and Its Impact on Crop Production

: Drought research is an important aspect of drought disaster mitigation and adaptation. For this purpose, we used the Standardized Precipitation Evapotranspiration Index (SPEI) to investigate the spatial-temporal pattern of drought and its impact on crop production. Using monthly precipitation (Precip) and temperature (Temp) data from 1986–2015 for 39 weather stations, the drought index was obtained for the time scale of 3, 6, and 12 months. The Mann–Kendall test was used to determine trends and rates of change. Precip and Temp anomalies were investigated using the regression analysis and compared with the drought index. The link between drought with large-scale atmospheric circulation anomalies using the Pearson correlation coefﬁcient (R) was explored. Results showed a non-uniform spatial pattern of dryness and wetness which varied across Myanmar agro-ecological zones and under different time scales. Generally, results showed an increasing trend for the SPEI in the three-time scales, signifying a high tendency of decreased drought from 1986–2015. The ﬂuctuations in dryness/wetness might linked to reduction crop production between 1986–1999 and 2005, 2008, 2010, 2013 cropping years. Results show relationship between main crops production and climate (teleconnection) factors. However, the low correlation values (i.e., <0.49) indicate the extent of the relationship within the natural variability. However, readers are urged to interpret this result cautiously as reductions in crop production may also be affected by other factors. We have demonstrated droughts evolution and trends using weather stations, thus providing useful information to aid policymakers in developing spatially relevant climate change adaptation and mitigation management plans for Myanmar. We examined how ENSO and IOD contributes to the drought evolution on correlation analysis performed at annual and seasonal scales. The study found that the ENSO results were positive and signiﬁcant mainly in the deltaic zone along the coastal region annual and August-October (ASO). While, the interior areas (i.e., the western hilly and central dry eco-zones) showed lower but negative correlation values for annual and ASO. Similarly, the IOD are similar (in R-values) but with opposite signs of the signiﬁcance level. In particular, we observe that eco-zones with highly signiﬁcant values in ENSO showed low to moderate values for IOD in annual and different seasons.


Introduction
The Intergovernmental Panel on Climate Change (IPCC) [1] stipulates that the worsening of global warming effects will lead to more extreme events such as droughts, floods, and heat waves. Evidence shows that some of these extreme events have already occurred in water-limited regions of the world [1]. Therefore, this study fills the gap by contributing to the literature and providing insight into drought studies and its relationship with crop production in Myanmar. FAO [50] reports have mentioned Asia economy to be highly vulnerable to climate change. As agriculture sector still remains the largest employer in most Asia economies, and its sensitivity to varying hydro-climatic conditions (particularly in regions with no or weak agricultural irrigation systems) cannot be overlooked.
In this study, we employ the SPEI developed by Vicente-Serrano et al. [22] based on the added value of representing the multi-scale characteristic of SPI and the sensitivity of PDSI to changes in evaporative demand. In addition, we estimate the link between weather and crop production for the three crops with the largest production value in the Myanmar: rice, wheat, and corn.
Myanmar is one of the most vulnerable countries in the region due to her diverse agro-climatic conditions [49]. Past studies showed that Myanmar has experienced both precipitation and temperature changes over the last decade [53][54][55][56]. These changes have had implication on the length of wet season and influenced farmers decision to grow their crops. For example, few studies found high precipitation and temperature variability to impact crop productivity [49,57]. However, the effect of climate change on Myanmar on the crop production is less explored.
We acknowledge that declines in staple crop production are caused by several factors including agricultural management practices, quality of seeds used, technological investments, fluctuation in crop prices, research and development [56], etc. Nonetheless, this study focuses on climate-driven crop production variability in Myanmar s major crops. The study outcomes may form a scientific basis to help design droughts mitigation and adaptation policies for different stakeholders in agriculture, water resources, and management, among others.
The present study aims at investigating the spatiotemporal variations, frequency of droughts and its possible connections with variations in SST and crop production anomalies using the SPEI at different time scales and seasons. In addition, it examines the difference between the SPEI based on the average of climatic parameters and the average of SPEI of all locations.
The remainder of this paper is organized as follows; provide a brief description of the study area and a summary of the methods and approaches used in Section 2. Section 3 provides the study results and discussion of findings. In Section 4. we provide the summary and conclusions of the study.

Study Area
Myanmar is situated in Southeast Asia. A variety of topography such as flat (central), hilly (north, northwest, and east), and coastal (west and south) covers the country. It is bounded by water (Bay of Bengal and the Andaman Sea) in the western part of the country with latitude 9 • 32 N-28 • [58]; (c) major rivers [58]; and (d) eco-physiographic zones [59].
Myanmar is characterized by precipitation a high spatio-temporal variability [53] as a result of its varying topography ( Figure 1a) and vegetation cover (Figure 1b). The vegetation cover interacts with the atmosphere in complex ways to affect weather and climate [60]. Based on land use and land cover (LULC) from DIVA-GIS [58], seven dominant classes of LULC were found over Myanmar, namely forest, open forest, degraded land, scrubland, agriculture land, mangroves, and water ( Figure 1b).
Myanmar has large rivers that crosses the country (Figure 1c). The figure presents the streamflow of four major rivers across Myanmar. The largest river is the Ayeyarwady (nearly 2,170 km in length), which cuts through the country′s fertile lands, followed by the Thanlwin, Sittaung and Chindwin Rivers. Those rivers define Myanmar′s agriculture sector contribution to Gross Domestic Product [50]. This explains the important role of agriculture in the region. More than 65% of the population live in rural areas, and employed under the agricultural sector [51]. The agro-ecological zone map of Myanmar is shown in Figure 1d.
Myanmar is characterized by precipitation a high spatio-temporal variability [53] as a result of its varying topography ( Figure 1a) and vegetation cover (Figure 1b). The vegetation cover interacts with the atmosphere in complex ways to affect weather and climate [60]. Based on land use and land cover (LULC) from DIVA-GIS [58], seven dominant classes of LULC were found over Myanmar, namely forest, open forest, degraded land, scrubland, agriculture land, mangroves, and water ( Figure 1b).
Myanmar has large rivers that crosses the country (Figure 1c). The figure presents the streamflow of four major rivers across Myanmar. The largest river is the Ayeyarwady (nearly 2,170 km in length), which cuts through the country s fertile lands, followed by the Thanlwin, Sittaung and Chindwin Rivers. Those rivers define Myanmar s agriculture sector contribution to Gross Domestic Product [50]. This explains the important role of agriculture in the region. More than 65% of the population live in rural areas, and employed under the agricultural sector [51]. The agro-ecological zone map of Myanmar is shown in Figure 1d.
By contrast, temperature variability is low, both spatially and temporally. The highest temperature (34.3 • C) is recorded in April till June, and the lowest (14.4 • C) is recorded in December till February [54]. The precipitation variability and frequent water scarcity unfavorably impact Myanmar s food security and socio-economic development agenda [61,62]. Moreover, atmospheric circulation drivers have been noted to influence Myanmar s and its surrounding inter-annual precipitation variability. For example, the El Niño Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD) events are linked to anomalous precipitation in the South East Asia region, leading to extreme floods and droughts, respectively.

In Situ Observation Data
We used the monthly mean precipitation and monthly mean temperature datasets of 39 meteorological stations spanning from 1986-2015. The data were collected from the Department of Meteorology and Hydrology, Myanmar. The location and specifics of each meteorological station are presented in Figure 1a and Table A1, respectively. The stations selection was selected based on their long-term temporal coverage, data homogeneity and data record's completeness. The SST indices such as El Niño-3.4 index derived from the Climatic Prediction Center (CPC) s, the Extended Reconstructed Sea Surface Temperature version 5 (ERSST v5) data [63]. The Dipole Mode Index (DMI) was derived from the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST1.1) [64] were used in this study.
To perform impact of climate change on crop production, the annual crop data series for Myanmar 1985-2015 were downloaded from the Food and Agriculture Organization Statistical Databases (FAOSTAT) website [65] (accessed June 6 2021) to analyze crop production anomalies. For crop production data, we selected three crops for further analysis; namely, corn, wheat and rice. These crops have sufficiently large acreage across Myanmar and uniform distribution over the arable land and thus are comparable during the study period. Figure S2 presents total crop production for the three major crops with their corresponding total hectares grown from 1984-2018.

Precipitation and Temperature Climatology
This study calculates monthly accumulated precipitation and temperature from daily data 1986-2015. Seasonal and annual values are then calculated from the monthly data.

Calculation of SPEI
To estimate PET, the Penman-Monteith and Thornthwaite models are the two most widely used models. The former provides a better estimate due to its more comprehensive physics [4,14] but it is data intensive and requires more input data that are not usually readily available. The latter is simple as it requires mean monthly temperatures.
The Thornthwaite method (Equation (1); [67] was used as it can capture the main impact of increased temperatures on water demand. The calculation method of Thornthwaite is as follows; Potential evapotranspiration (PET) was calculated using Equation (1).
where T is the monthly mean temperature ( • C); I is a heat index calculated for the whole year; m is a coefficient depending on I, such that m = 6.75 × 10 − 7I 3 − 7.71 × 10 − 5 I 2 2 + 1.79 × 10 − 2I 2 + 0.492; and K is a correction coefficient defined based on the latitude and month. The difference between the precipitation (precip) and PET is computed using Equation (2) [22].
D is the difference between the precipitation Precip and PET for the month i and PRECIP is the precipitation. The difference D k i,j in a given month j of year i depends on the Agronomy 2021, 11, 1691 6 of 28 timescale k. For example, the accumulated difference for month i in a particular year i with a 12-month timescale is computed using Equation (3).
where D I,I (mm) is the PRECIP −PET difference in the first month of the year i. The water balance is normalized based on a log-logistic probability distribution to obtain the SPEI time series. The probability density function of a variable with a log-logistic distribution is expressed as Equation (4): where α, β, and γ are the scale, shape, and origin parameters, respectively, for D values in the range γ < D < ∞. Thus, the probability distribution function of series D f(x) is given by the following Equation (5).
The SPEI was estimated in an R statistical software [68]. Drought was analyzed with SPEI at 3, 6, 12-month timescales [19]. To better understand the extent and intensity of drought, the drought frequency (F) value was used. The F-value was calculated using the following Equation (6).
where, n is the number of drought years at the station, N is the total number of years for the precipitation data, and i represents the meteorological station.

Standardized Anomaly Estimation
We estimated the standardized anomaly for precipitation, temperature and the three different crops production in Myanmar based on Equation (A6) in Appendix A.3.

Statistical Analysis
Mann-Kendall test [69] was used to determining trends and rates of change. The formula and estimation procedure are presented in Equations (A1)-(A4) in Appendix A.1. In addition, we used the least square method Equation (A5) to compute the linear trend in precipitation and temperature anomalies. We also performed a correlation analysis Equation (A6); Appendix A.2) to describe the relationship between SPEI and remote drivers (i.e., ENSO and IOD).
In addition, performed a correlation analysis to assess the link between annual SPEI and crop production across multiple years. We also analyzed using seasonal SPEI to investigate its impact on crop production since Myanmar has two crop growing seasons. We selected May-June-July (MJJ), August-September-October (ASO), and MJJASO as seasons to study its effect on crop production.
Statistical analysis was performed using R and Python softwares [68]. Many hydrometeorological studies have adopted this kind of approach [70][71][72][73].  Figure 2 presents the long-term spatial distribution (1986-2015) annual mean of Temp and Precip. Temp has a bi-modal distribution with warmer temperatures on central, Rakhine coastal, Yangon deltaic, Ayeyarwady deltaic and southern Myanmar coastal and cooler temperatures on the north hilly and eastern hilly (Figure 2b). The results of Temp explain the country received the coldest place on north hilly and eastern hilly while hottest place on central, Rakhine coastal, Yangon deltaic, Ayeyarwady deltaic and southern Myanmar coastal (Figure 2b).  Figure 2 presents the long-term spatial distribution (1986-2015) annual mean of Temp and Precip. Temp has a bi-modal distribution with warmer temperatures on central, Rakhine coastal, Yangon deltaic, Ayeyarwady deltaic and southern Myanmar coastal and cooler temperatures on the north hilly and eastern hilly (Figure 2b). The results of Temp explain the country received the coldest place on north hilly and eastern hilly while hottest place on central, Rakhine coastal, Yangon deltaic, Ayeyarwady deltaic and southern Myanmar coastal (Figure 2b). On the other hand, Precip and Temp follows the topography. The hilly region received less Temp and coastal and central dry zone received high Temp. The Rakhine coastal and southern Myanmar coastal received highest Precip due to the modulated by summer monsoon season over Southeast Asia.

Climatology and Linear Trend of Temperature and Precipitation
Also, the PET results indicate that the effect of PET is an important hydro-climatic variable in the Myanmar ( Figure S1). The PET values are linked with the Precip and Temp distribution in the region. The spatial pattern of PET ( Figure S1a) results shows the north hilly and eastern hilly receive low PET values. However, central, Rakhine coastal, Yangon deltaic, Ayeyarwady deltaic and southern Myanmar coastal show high PET.  On the other hand, Precip and Temp follows the topography. The hilly region received less Temp and coastal and central dry zone received high Temp. The Rakhine coastal and southern Myanmar coastal received highest Precip due to the modulated by summer monsoon season over Southeast Asia.
Also, the PET results indicate that the effect of PET is an important hydro-climatic variable in the Myanmar ( Figure S1). The PET values are linked with the Precip and Temp distribution in the region. The spatial pattern of PET ( Figure S1a) results shows the north hilly and eastern hilly receive low PET values. However, central, Rakhine coastal, Yangon deltaic, Ayeyarwady deltaic and southern Myanmar coastal show high PET. Figure 3 shows the annual variation in Temp and Precip over Myanmar from 1986-2015. We observed that the Temp pattern is bi-modal, with the highest peak in April (summer season) and the second in October (late monsoon period). Moreover, a unimodal precipitation pattern corresponding to the summer monsoon season from May-October is observed. The highest Precip peak is found in June-July-August. Overall, we found a distinctive characteristic in spatial Temp and Precip patterns ( Figure 3). When the Precip is high as oppose Temp is low in the region. It is modulated by Indian monsoon season over Myanmar. In addition, the estimated values in temporal PET patterns ( Figure S1b) are largely associated with the Precip, Temp and other hydro-climate and teleconnections variabilities in the region.
Agronomy 2021, 11, x FOR PEER REVIEW 8 of 28 shows the annual variation in Temp and Precip over Myanmar from 1986-2015. We observed that the Temp pattern is bi-modal, with the highest peak in April (summer season) and the second in October (late monsoon period). Moreover, a unimodal precipitation pattern corresponding to the summer monsoon season from May-October is observed. The highest Precip peak is found in June-July-August. Overall, we found a distinctive characteristic in spatial Temp and Precip patterns ( Figure 3). When the Precip is high as oppose Temp is low in the region. It is modulated by Indian monsoon season over Myanmar. In addition, the estimated values in temporal PET patterns ( Figure S1b) are largely associated with the Precip, Temp and other hydro-climate and teleconnections variabilities in the region. In addition, we plotted the inter-annual variation in precipitation and temperature standardized anomalies from 1986 to 2015 ( Figure 4). The result showed both temperature and precipitation exhibited an increasing trend. In addition, we plotted the inter-annual variation in precipitation and temperature standardized anomalies from 1986 to 2015 ( Figure 4). The result showed both temperature and precipitation exhibited an increasing trend. Figure 5 describes the area average of SPEI (SAA) and SPEI based on area average of climatic parameters (SBA) over the period of study . The SAA was obtained by averaging SPEI of each location, while SBA is directly obtained after averaging the monthly precipitation and temperature before calculating the SPEI over Myanmar. The SAA was standardized to unit before calculating its interdependence with SBA. Significant and positive correlation coefficients are obtained between SAA and SBA for the three-time scales, i.e., 3-, 6-and 12-months. This means that area average of climate variables is useful to characterize the drought condition over Myanmar, which is in line with Barren et al., [74].

Temporal Variations in Wetting and Drying Trends
Overall, the temporal patterns for SPEI-3, SPEI-6 and SPEI-12, and the patterns for years with dryness and wetness are similar however differences in their relative magnitudes can be noted. It is found that short-term climate conditions (SPEI-3 and SPEI-6) in Figure 5a,b are closer to each other than that of long-term climatic conditions (SPEI-12).
From Figure 5a, it is observed that the SPEI-3 and SPEI-6 and SPEI-12 showed 1998 as the extreme drought year. Meanwhile, the magnitude of drought intensity are higher as the time scale increases.   Figure 5 describes the area average of SPEI (SAA) and SPEI based on area average of climatic parameters (SBA) over the period of study . The SAA was obtained by averaging SPEI of each location, while SBA is directly obtained after averaging the monthly precipitation and temperature before calculating the SPEI over Myanmar. The SAA was standardized to unit before calculating its interdependence with SBA. Significant and positive correlation coefficients are obtained between SAA and SBA for the threetime scales, i.e., 3-, 6-and 12-months. This means that area average of climate variables is useful to characterize the drought condition over Myanmar, which is in line with Barren et al., [74].  Table 1 shows the classification of SPEI s (in terms of intensities) over Myanmar from 1986 to 2015. The SPEI-3 showed the highest occurrence of droughts, followed by SPEI-6. SPEI-12 showed the lowest percentage of extreme droughts. Between 1986 and 2015, it is found that SPEI-3 show dominant severe droughts over very wet conditions. Meanwhile, the SPEI-6 show a dominant very wet conditions over severe droughts.

Temporal Variations in Wetting and Drying Trends
Concerning the SPEI-12, a high occurrence of extreme drought is observed, while a low occurrence of extreme wet condition is evident.  Overall, the temporal patterns for SPEI-3, SPEI-6 and SPEI-12, and the patterns for years with dryness and wetness are similar however differences in their relative magnitudes can be noted. It is found that short-term climate conditions (SPEI-3 and SPEI-6) in Figure 5a,b are closer to each other than that of long-term climatic conditions (SPEI-12).
From Figure 5a, it is observed that the SPEI-3 and SPEI-6 and SPEI-12 showed 1998 as the extreme drought year. Meanwhile, the magnitude of drought intensity are higher as the time scale increases. Table 1 shows the classification of SPEI′s (in terms of intensities) over Myanmar from 1986 to 2015. The SPEI-3 showed the highest occurrence of droughts, followed by SPEI-6.

Spatial Variations in Wetting and Drying Trends
The monotic trend analysis of SPEI at 3, 6, 12-months was examined using the Mann-Kendall test at 5% significance level [69]. Figure 6 shows the precipitation trends from 1986 to 2015. These results help us understand the spatial distributions of individual station linear trends for different SPEI during the study period.

Spatial Variations in Wetting and Drying Trends
The monotic trend analysis of SPEI at 3, 6, 12-months was examined using the Mann-Kendall test at 5% significance level [69]. Figure 6 shows the precipitation trends from 1986 to 2015. These results help us understand the spatial distributions of individual station linear trends for different SPEI during the study period. Results of SPEI-3 show increasing trends mainly in the western hilly, deltaic (along the coastal region), and central and east parts (Figure 6a). We observed that few stations located in the central dry area show increasing and decreasing trends. On the contrary, stations located in the northern hilly zone show a decreasing trend (Figure 6a). Results of SPEI-3 show increasing trends mainly in the western hilly, deltaic (along the coastal region), and central and east parts (Figure 6a). We observed that few stations located in the central dry area show increasing and decreasing trends. On the contrary, stations located in the northern hilly zone show a decreasing trend (Figure 6a).
The spatial pattern of the SPEI-6 shows increasing trends at many stations in the deltaic zone (but not significant at both 5 and 10%, except for few stations). We observe larger increasing trends in the southern-most portion of the deltaic zone.
Notably, in the central hilly and western hilly, many stations showed a mixed spatial pattern of increasing and decreasing trends. The SPEI-6 results for the northern hilly area are similar to SPEI-3, with the stations showing decreasing trends. SPEI-12 denotes longterm droughts (Figure 6c). The SPEI-12 results are similar to SPEI-6 but with statistically significant results across Myanmar. Moreover, the highest increasing trends are shown for long-term droughts.
We quantitively evaluated the spatial variability of drought frequency across Myanmar in the different SPEI by computing the SPEI frequency distribution for the annual and seasonal basis to understand the regional wetting and drying conditions. We obtained the drought frequency at each station and based on a spatial analysis method, interpolated the results to obtain different drought frequencies at a country-level. Figures 7 and 8 presents the spatial distribution of drought frequency for the three different levels of drought for 1986-2015.
We observed that the western hilly areas and patches in the deltaic eco-zones along the coastal region showed the highest frequency of dry spells as compared with wet spells in the same period. The findings showed similar spatial distributions for different seasons but with varying values (Table 1). mar in the different SPEI by computing the SPEI frequency distribution for the annual and seasonal basis to understand the regional wetting and drying conditions. We obtained the drought frequency at each station and based on a spatial analysis method, interpolated the results to obtain different drought frequencies at a country-level. Figures 7 and 8 presents the spatial distribution of drought frequency for the three different levels of drought for 1986-2015.  Taking MJJ as an example, fewer dry spells than wet spells persisted with varying frequencies (Figures 7b and 8b). In ASO, the frequency of dry spells is found in nearly the entire country, with deltaic, western, and eastern hilly zones showing a more pronounced frequency of dry spells. In contrast, wet spells in ASO (Figure 8c) showed similar spatial patterns to Figure 7c (MJJ wet spells), but with significant frequency values. In the monsoon season, the frequency of dry spells was high (4-7 times) in most areas from deltaic, western, and eastern hilly zones compared with wet spells in both space and time.
frequencies (Figures 7b and 8b). In ASO, the frequency of dry spells is found in nearly the entire country, with deltaic, western, and eastern hilly zones showing a more pronounced frequency of dry spells. In contrast, wet spells in ASO (Figure 8c) showed similar spatial patterns to Figure 7c (MJJ wet spells), but with significant frequency values. In the monsoon season, the frequency of dry spells was high (4-7 times) in most areas from deltaic, western, and eastern hilly zones compared with wet spells in both space and time.

Correlation Analysis of SPEI and Its Key Mechanisms
To understand the underlying remote drivers contributing to the drought evolution across Myanmar, we examine the link between droughts evolution and some key atmospheric circulation drivers. Here, we focused on two drivers, namely, El Niño-Southern

Correlation Analysis of SPEI and Its Key Mechanisms
To understand the underlying remote drivers contributing to the drought evolution across Myanmar, we examine the link between droughts evolution and some key atmospheric circulation drivers. Here, we focused on two drivers, namely, El Niño-Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD), based on a recommendation from past studies [75], to understand the spatial and temporal dynamics.
To calculate the spatial pattern, correlation analysis was performed at annual and seasonal scales. Figure 9 and Figure S1 present correlation and their significance between droughts, ENSO and IOD, respectively.
To calculate the spatial pattern, correlation analysis was performed at annual and seasonal scales. Figure 9 and Figure S1 present correlation and their significance between droughts, ENSO and IOD, respectively.
In Figure 9a, we observe the correlation coefficients between SPEI the annual values and Niño-3.4 index. The correlation results were positive and significant (i.e., ≥ 0.58, p < 0.05) mainly in the deltaic zone along the coastal region.   The yellow arrow denotes no correlation; purple denotes weak correlation; the light green denotes moderate while the deep green denotes high correlation, respectively. The blue and red dots represent wet and dry, respectively.
In Figure 9a, we observe the correlation coefficients between SPEI the annual values and Niño-3.4 index. The correlation results were positive and significant (i.e., ≥ 0.58, p < 0.05) mainly in the deltaic zone along the coastal region. Figure 9b-d show station correlation for individual seasons. In general, the spatial patterns for annual and August-October (ASO) show similar patterns but differences in levels of significance. However, we observed that the spatial patterns for MJJ and the entire monsoon season (MJJASO) were similar in both magnitude of correlation coefficients and significance levels.
The western hilly and central dry eco-zones showed lower but negative correlation values for annual (i.e., ≥0.25, p < 0.05) followed by ASO (≤0.42, p < 0.05) and mixed results in the northern and eastern hilly eco-zones. The mixed results in those regions range from high positive (and highly significant) to moderately negative (significant).
In Figure S1a-d, the correlation coefficients between annual ( Figure S1a) and seasonal SPEI ( Figure S1b-d) and IOD are similar (in R-values) but with opposite signs of the significance level. In particular, we observe that eco-zones with highly significant values in ENSO showed low to moderate values for IOD in annual and different seasons.
Secondly, we performed correlation at each station with ENSO and IOD on an annual scale to understand the temporal dynamics. Table 2 shows the correlation (R) between SPEI and ENSO and IOD. The results show that the SPEI influence was mainly explained by ENSO and IOD episodes. We observe that an overall weak correlation value (no significant level) concerning ENSO (≥0.15, p > 0.05) while in IOD results (≥0.46, p < 0.05), the correlation coefficient results moderate and statistically significant at 5% level. Generally, the SPEI influence was mainly explained by changes in IOD in the region since the correlation values were reported at 0.46 of SPEI variations ( Table 2).

Crop Productions and Its Influencing Factors
To describe the major agriculture characteristics of Myanmar, we present the temporal distribution of the three major crops grown in Myanmar. The crop annual statistics data were downloaded from the FAOSTAT database [65]. The annual distributions of the major three crops grown in Myanmar from 1986 to 2015 are presented in Figure 10. Figure 10a presents the production and total area statistics of rice in Myanmar. Generally, from 1984 to 2018, rice shown a steady increase in production ( Figure 10, red line graph) and a corresponding increase in rice area grown ( Figure 10, bar graph). We observe steady rice production (area grown) move from 7075 × 10 3 MT (4603 HA) in 1986 to 13,200 × 10 3 MT (7040 HA) in 2015.
The result in Figure 10 makes rice as the first extensive crop grown followed by corn production in Myanmar (Figure 10a). This result is not surprising as rice production the most geographically ubiquitous crop in Myanmar (as well as south and southeast Asia). In addition, the Figure 10b result makes corn is the second major crop grown in Myanmar from 1986 to 2015 consistent with FAO reports [76].
In addition, to understand possible climate changes implications on crop production, we conducted two separate analyses; the composite and correlational analysis. Reasons for selecting these two analyses is well articulated in [4,77]. Here, the study explored how to combined historical observations of crop production and weather and climate variables to explains crop production. The period of our analysis was based on annual timescales from 1986 to 2015 consistent with the period of weather station acquired.
The result in Figure 10 makes rice as the first extensive crop grown followed by corn production in Myanmar (Figure 10a). This result is not surprising as rice production the most geographically ubiquitous crop in Myanmar (as well as south and southeast Asia).  Figure 11 shows the composite analysis of rice production anomalies with multiple influencing factors (i.e., SPEI, precipitation anomalies, temperature anomalies, teleconnections).
Generally, we observed a significant decrease in crop production from 1986 to 1999 for rice (Figure 11a), corn ( Figure 12a) and Wheat ( Figure S3a) as compared to 2000-2015. However, we observed that two contrasting periods stood out (i.e., 1986-1999 and 2002-2015) consistent to relatively, the period of dryness (wetness) from 1986-1999 (2000-2015) in Figures 4a and 5c. from 1986 to 2015 consistent with the period of weather station acquired.
For ease of comparison, the study made a composite of annual SPEI (in Figure 5c), precipitation (Figure 4a), temperature (Figure 4b) anomalies to present as shown in Figure  10b, Figure 11b and Figure S3b while the remaining figures were new information. Figure 11 shows the composite analysis of rice production anomalies with multiple influencing factors (i.e., SPEI, precipitation anomalies, temperature anomalies, teleconnections). From both Figures 4a and 5c we observed similar dryness periods (1986-1989, 1992/3, 1995/6, 1998, 2003-2005) varying from mild to moderate drought (with values <−1.5) and relatively wetness periods in the remaining periods. The most severe drought event on record occurred from 1998,2005,2008,2010,2013 in Figure 5 corroborates with the period of crop production reduction (Figure 11, Figure 12 and Figure S3). Also, Figure 11b,c, Figure 12b,c and Figure S3b,c) showed the composite analysis of El Niño episodes (2001,2002,2003) and La Niña episodes (2010-2013 which occurred over three conservative ENSO years within the period 1986 to 2015, which has implications on atmospheric water supply (i.e., precipitation) and heat stress (temperature rise) over the regions. This result is consistent with the previous literature [78,79]. Similar results across the globe are reported in the IPCC AR4&5 reports [1,80] with a corresponding threat to different sectors of the society (e.g., agriculture, economic and social) temperature, red-SPEI) (c) DMI anomalies, (c) Nino 3.4 SST anomalies during the same time period.
Generally, we observed a significant decrease in crop production from 1986 to 1999 for rice (Figure 11a), corn (Figure 12a) and Wheat ( Figure S3a) as compared to 2000-2015. However, we observed that two contrasting periods stood out (i.e., 1986-1999 and 2002-2015) consistent to relatively, the period of dryness (wetness) from 1986-1999 (2000-2015) in Figures 4a and 5c. From both Figures 4a and 5c we observed similar dryness periods (1986-1989, 1992/3, 1995/6, 1998, 2003-2005) varying from mild to moderate drought (with values < −1.5) and In addition, crops growth is sensitive to different aspects of the climate and under the current warming climate, this study explored the relative importance of different climatic factors on the three-crop production across Myanmar. We conducted correlational to helps us to understand which specific crop types were significantly exposed to dryness and wetness conditions following similar analysis in [81]. The main difference in our study analysis from [71], is the including of multiple datasets (i.e., long term (annual) and short (i.e., seasonal (MJJ, ASO and MJJASO)).
The selection of the seasons was done following the three cropping seasons (summer, monsoon, and winter) of Myanmar. The summer season goes February-May, while monsoon is normally June-September, and winter is October-January [49,82]. In Table 3, we presented the correlation results of crop anomalies with the multiple datasets. Table 3. Correlation coefficients correlation coefficients between crop production anomalies (CPA) and climate and teleconnection variables from 1986-2015. From Table 3, the results showed that the relative importance of each of the climate variables vary in their response to crop production anomalies. We observed that the correlation values ranging from −0.009 to 0.460 (Table 3). Although the R values were relatively low this does not mean that climate variables in Myanmar do not have a trend, but just that the trend stays within the limits of natural variability during the time range of our analysis. This assertation is consistent studies conducted at global scale [83] and regional scales [46,84].
Generally, we found that SPEI (annual) the accounted for 0.12 of the combined effect on the three-crop production variability. This suggest that over the period, the variations of these three-crop production are explained by changes in SPEI in the region. However, the individual coefficient values varied considerably with Corn (R = 0.280) showed higher variability, followed by Rice (R = 0.02) and Wheat (0.02).
On the other hand, when we examined the shorter time scales drought indices (MJJ, ASO), we observed that the average value increased in MJJ (R = 0.147) and ASO (R = 0.160) for the three crops but lower for MJJASO (R = 0.080) suggesting that the crops are most slightly correlated in MJJ and ASO than on longer time scales. This result corroborates with previous studies that suggest that shorter time scale is appropriate for agricultural drought monitoring [85].
Also, we found among the three crops, corn production anomalies increase in MJJ (R = 0.320) and ASO (R = 0.390) but lower in MJJASO (R = 0.170). The results for MJJ and ASO for corn production are not surprising these months are the most critical cropping months (i.e., the early monsoon season) which correspond to the corn planting to early reproductive stage (tasseling/silking) in Myanmar cropping calendar. The result is consistent with previous studies in the region [82] as the phenology of these crops explains why their anomalies correlate most closely with MJJ than the rest. This result is consistent with findings of [69,71,72].
On the contrary, wheat showed that positive correlation in MJJ (R = 0.160) and negative R = −0.03 in ASO. The ASO is late monsoon season in Myanmar and correspond to harvesting time for wheat. Rice result showed that negative correlation in MJJ (R = −0.04) and positive R = 0.12 in ASO.
We observe precipitation, and temperature anomalies show the highest correlation with corn production anomalies, followed by Wheat and Rice ( Table 3).
The PRECIP positive relationship was considerably weaker (R = 0.046) than in TEMP (R ≤ 0.110) ( Table 3) in the region suggesting that TEMP is a limiting driver of crop production anomalies than PRECIP in the region.
Crops are known to be impacted by precipitation and temperature variability [86]. Table 3 shows corn production anomalies showed the highest correlation values (in both PRECIP and TEMP) than Wheat and Rice. Note that rice shows a negative but weaker correlation of with TEMP (Table 3, R = −0.09). Similar results were shown for corn for SPEI at both annual and seasonal scales (Table 3). This result is not surprising as corn is sensitive to water stress (i.e., both water supply and demand, respectively) at five main phenological stages of corn growth, although their sensitivity varied at each stage [87].
In Table 3, the results of wheat showed responds to TEMP (R = 0.34) than PRECIP (R = 0.09). These results suggest respond to heat stress (i.e., TEMP) than PRECIP. This result finding is in line with Semenov and Shewry [88]. In addition, the wheat relationship with TEMP is not surprising as this result is consistent with a study by Lobell et al. [46] in India.
Generally, crops responds to TEMP, when coincides with flowering, impacting the ageing of leaves and accelerates crops towards maturity consistent with [46,79]. Schelenker and Roberts [84], found high TEMP influences crops yields in the US. TEMP may be a more significant driver of crop production than precipitation. However, more detailed analysis is needed to understand this mechanism at field scale. Further studies should use longer time period and/or grided datasets to confirm this analysis ENSO and IOD play key role in the crop production cycle. Results of CPA with ENSO and IOD is shown in Table 3. ENSO and IOD had a significant positive correlation coefficient (R) with corn (R = 0.203) but negative with wheat (R = −0.133) and rice (R = −0.004), respectively. This result agrees well with that ENSO and IOD cycle impacts crop via controls on precipitation. This result corroborates with previous studies that suggest ENSO and IOD cycles impacts crop productions [89,90]. We acknowledge that quantifying and comparing the relationship between drought losses across time is challenging as crop productions is not controlled only by weather and climate factors but by other factors such as scientific and technological advances (e.g., improvements in plant genetics, fertilizer, pesticides, and irrigation facilities). In addition, the spatial sampling of the meteorological stations could explain the low correlation values. For example, the meteorological stations at a given point is not as representative of the entire area. The spatial sampling of the localized rain cells may occasionally poorly capture precipitation event to influence the correlation results. This assertation is consistent with [91]. We recommend future studies may consider using grided datasets to examine this phenomenon. However, overall, this study has demonstrated that the fluctuations of crop production may be explained by weather and climate factors.

Summary and Conclusions
The long-term spatial distribution (1986-2015) annual mean of Temp and Precip was analyzed. Overall, the Precip and Temp follows the topography. The hilly region received less Temp and coastal and central dry zone received high Temp. The Rakhine coastal and southern Myanmar coastal received highest Precip due to the modulated by summer monsoon season over Southeast Asia. The annual variation in Temp and Precip over Myanmar showed that the Temp pattern is bi-modal, with the highest peak in April (summer season) and the second in October (late monsoon period). Moreover, a unimodal precipitation pattern corresponding to the summer monsoon season from May-October is observed. The highest Precip peak is found in June-July-August. A distinctive characteristic in spatial Temp and Precip patterns (Figure 3) showed that the high Precip as oppose low Temp in the region is modulated by Indian monsoon season over Myanmar.
Also, the inter-annual variation in Precip and Temp standardized anomalies presented an opposite variation between the two variables ( Figure 4). The high (low) precipitation value coincides with low (high) temperature value, however, overall, both variables exhibited an increasing trend. The analysis of SPEI over Myanmar was presented for multiple time scales of 3,6, and 12 months, respectively ( Figure 5). Overall, the temporal patterns for SPEI-3, SPEI-6 and SPEI-12, and the patterns for years with dryness and wetness are similar however differences in their relative magnitudes can be noted. It is found that short-term climate conditions (SPEI-3 and SPEI-6) in Figure 5a,b are closer to each other than that of long-term climatic conditions (SPEI-12) with 1998 recorded as an extreme drought year in the SPEI-3 and SPEI-6 and SPEI-12. On the other hand, the SPEI-3 showed the highest occurrence of droughts, followed by SPEI-6. SPEI-12 showed the lowest percentage of extreme droughts in terms of drought intensity (Table 1). SPEI-3 show dominant severe droughts over very wet conditions between 1986 and 2015.
We examined how ENSO and IOD contributes to the drought evolution on correlation analysis performed at annual and seasonal scales. The study found that the ENSO results were positive and significant mainly in the deltaic zone along the coastal region annual and August-October (ASO). While, the interior areas (i.e., the western hilly and central dry eco-zones) showed lower but negative correlation values for annual and ASO. Similarly, the IOD are similar (in R-values) but with opposite signs of the significance level. In particular, we observe that eco-zones with highly significant values in ENSO showed low to moderate values for IOD in annual and different seasons.
Within the agricultural sector, droughts reduce soil-water availability, thus, affecting crop failures and pasture losses. If drought occurrences last longer then it may severely impact growth of crops. There is evidence to show that these severe droughts have impacted crop yield thus, threatening regional food security in regions with rain-fed agricultural system. The potential influence of climate variability on crop production would be its sensitivity mainly to rainfall and temperature variability.
The study used both composite and correlational analysis to understand possible climate changes implications on crop production. The composite analysis of crop production anomalies analyzed with multiple influencing factors (i.e., SPEI, precipitation anomalies, temperature anomalies, teleconnections). A significant decrease in crop production from 1986 to 1999 was observed for all crops compared to 2000-2015 consistent with the relative dryness and wetness conditions over the period. The two contrasting periods stood out (i.e., 1986-1999 and 2002-2015) consistent to the period of dryness (wetness) from 1986-1999 (2000-2015). Note that these periods were characterized with dryness periods (1986-1989, 1992/3, 1995/6, 1998, 2003-2005) varying from mild to moderate drought (with values <−1.5) and relatively wetness periods in the remaining periods. Similarly, the composite analysis of El Nino episodes (2001,2002,2003) and La Niña episodes (2010-2013 which occurred over three conservative ENSO years within the period 1986 to 2015, which has implications on atmospheric water supply (i.e., precipitation) and heat stress (temperature rise) over the regions.
The correlational helps us to understand which specific crop types were significantly exposed to dryness and wetness conditions. Crop growth is sensitive to different aspects of the climate and under the current warming climate. Based on the three cropping seasons (summer, monsoon, and winter) of Myanmar, the correlation results of crop anomalies with the multiple datasets showed interesting results. The results showed that the relative importance of each of the climate variables vary in their response to crop production anomalies. The correlation values ranging from −0.009 to 0.460. Although the R values were relatively low this does not mean that climate variables in Myanmar do not have a trend, but just that the trend stays within the limits of natural variability during the time range of our analysis. This assertation is consistent studies conducted at global scale and regional scales.
Overall, the variations of these three-crop production are explained by changes in the multiple datasets used in the region. Here, our analysis focused on how the climate change signal on crop production anomalies from 1985 to 2015 to understand with sufficient detail the processes involved in potential crop production changes. We urge readers to interpret this result findings caution as crop production declines are not driven by a single event, but rather result from a confluence of other factors including farm management practices adopted, seed quality used, technological investments inputs, crop price fluctuations, research and development, etc. We have identified the spatiotemporal variation of drought and its impacts on agriculture. This results analysis could provide policy makers and stakeholders with scientific information regarding which agricultural areas are most vulnerable and sensitive to drought. We recommend further studies could build upon our analysis. Future studies on trend analysis should consider non-linear models that are predictive and also include physiological detail necessary to explain the processes behind the crop production-climate (teleconnection) associations. In addition, further study should study the future changes in climate and its impact on food production over the region. Particularly increasing resource constraints within the context of climate change, with decreasing water and land resources would impact food production.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11091691/s1, Figure S1. Spatial distribution of (a) PET (mm) climatology and (b) Annual cycle of PET (mm) over Myanmar from averaged from 1986 to 2015. Figure S2 Same as Figure 9, but for correlation between IOD and SPEI. Figure

Acknowledgments:
We are grateful to the two anonymous reviewers for their time and support in reshaping the manuscript. Their comments helped us to improve the final version of the manuscript. Furthermore, this research was encouraged by Binjiang College, Nanjing University of Information Science and Technology, Wuxi City, Jiangsu Province, China. Special appreciation goes to the Department of Meteorology and Hydrology, Myanmar, to provide the datasets used in the study.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A.1 Calculation Trends Analysis
Trends in drought variations were calculated using the Mann-Kendall tau-b nonparametric technique. This approach is widely adopted for hydro-meteorological time series [92]. The tests were calculated based on Equations (6) and (A1)-(A3).
where S is the rating score (called the Mann-Kendall sum); x is the data value; I and j counters; n represents the number of data values in the series; sgn x i − x j is a function. Positive and negative values of S indicate increasing and decreasing trends, respectively. The variance is Var (S) = n(n − 1)(2n + 5) 18 (A2) S is standardized (as shown in Equation (A2)) by subtracting its expectation (zero) divided by its standard deviation (σ S ). , i f S < 0 (A3) |Z| > Z ∝/2 signifies that the time series data show a significant trend. ∝ is the significance level. The present study sets the significance level to 0.05, corresponding to Z ∝/2 = 1.96. Thus, when the time series data produce |Z| > 1.96, there is a significant increase or decrease trend.

Appendix A.2 Linear Regression Model and Correlational Analysis
A linear regression Equation (A5) was employed to calculate the linear trend in monthly and seasonal monsoon precipitation.
where, Y represents the dependent variable (i.e., precipitation, temperature, crop productions trends, respectively); β o and β 1 are the coefficients. Time (t) is the predictor. A Pearson correlation coefficient (R 2 ) Equation (A5) was used to examine the effects of either ENSO and IOD on droughts variations.
where n represents the number of months and Ai and Bi represent the monthly SPEI and the teleconnections (i.e., ENSO and IOD) and crops (corn, wheat and rice) data at time i, respectively. However, before computing the correlation function between crop productions and droughts, we detrended the crop production data using a specific predetermined function model following [93][94][95].

Appendix A.3 Standardised Anomaly
The standardized rainfall values were calculated for all the years from the long-term mean, yearly mean, and the standard deviation using Equation (A6): where Z represents the standardized departure, X is the variable at a certain period (i), X is the long-term mean value, and S d is the standard deviation from the mean. The Z value provides immediate information about the significance of a specific deviation from the mean [96].