Trends in Long-Term Drought Changes in the Mekong River Delta of Vietnam

In recent years, short droughts in the dry season have occurred more frequently and caused serious damages to agriculture and human living in the Mekong River Delta of Vietnam (MRD). The paper attempts to quantify the trends of drought changes in the dry seasons from 2001 to 2015 in the region, using daily MODIS MOD09GQ and MOD11A1 data products. Here, we exploit the Temperature Vegetation Dryness Index (TVDI) to assess levels of droughts. For each image-acquisition time, the TVDI image is computed, based on the Normalized Difference Vegetation Index (NDVI), derived from red and near infrared reflectance data, and the Land Surface Temperature (LST), derived from thermal infrared data. Subsequently, a spatiotemporal pattern of drought changes is estimated, based on mean TVDI values of the dry seasons during the observed period, by a linear regression. As a result, the state of drought in the dry seasons in the MRD has mostly been at light and moderate levels, occupying approximately 62% and 34% of the total area. Several sub-areas in the center have an increased trend of drought change, occupying approximately 12.5% of the total area, because impervious surface areas increase, e.g., the obvious land use change, from forest land and land for cultivation for perennial trees being strongly converted to built-up land for residence and public transportation. Meanwhile, several sub-areas in the coastal regions have a negative trend of drought change because water and absorbent surface areas increase, e.g., most of land for cultivation for perennial trees has been converted to aquaculture land. These cases usually occur in and surrounding forest and wet land, also occupying approximately 12.5% of the total area.


Introduction
Drought is a natural phenomenon that seriously affects agricultural production and human living. Recently droughts occur in many countries or regions with significant impacts and increase in frequency, severity, and duration [1][2][3]. According to the National Oceanic and Atmospheric Administration (NOAA), global drought area reached its highest level in 2015 to 2016 during the last decades [2]. In the past, drought is driven by natural variation in seasonal or annual precipitation, so the risk of drought is expected to grow due to reduced precipitation and higher temperature [4]. In recent years, the severity of a drought event can be affected by human-induced climate change in combination with natural variations [4,5]. Additionally, increased heating from global warming is considered to set droughts occurring quicker and more intense [4][5][6]. Therefore, monitoring and assessing drought events is an actual challenge worldwide.

Materials and Methods
In this section, we firstly introduce the study area, the MRD. Secondly, we describe MODIS input data, including MOD09GQ and MOD021KM products. Then, we conduct a data processing to compute TVDI values. Finally, a trend analysis of TVDI changes in the dry seasons between 2001 and 2015 are estimated for exploring drought changes in the MRD.

Study Area
The MRD is located in the southwestern part of Vietnam and is the ending part of the Mekong River Basin, as shown in Figure 1. It borders to Cambodia in the north and northwest, Gulf of Thailand in the west and southwest, and the South China Sea in the east and southeast. It occupies the total area of 40,816 km 2 [43]. Except for a mountainous area in the north, consisting of seven discontinuous mountains with the highest 705 m peak above sea level and the slope of more than 25 degrees, belonging to the two Tinh Bien and Tri Ton districts of the An Giang Province, the region has rather flat terrain and is at low altitude above sea level.
The climate in the MRD is tropical, hot, and humid, and dominated by the Asian monsoons. According to General Statistics Office of Vietnam [43], its mean daily temperature is variable through the year from 24 • C to 29 • C, while in a year the highest temperature is about 36.3 • C and the lowest is in the order of 18.0 • C. The sunshine hours range from about 2000 to 2500 h per year, occurring much more frequently in the dry season than in the wet season. The MRD has an annual rainfall of more than 2000 mm, but heavy showers regularly occur during the wet season. For example, at the Ca Mau meteorological station from 2002 to 2015, the monthly averages of air temperature, sunshine hour, and precipitation are about 27.3 • C, 210 h, and 40 mm in the dry season, respectively, while those are about 27.8 • C, 155 h, and 300 mm in the wet season [43]. Additionally, a no or very low precipitation between January and March normally occurs over the MRD, also see Section 4.1 for an analysis of precipitation at the Can Tho and Dong Thap meteorological stations. Therefore, the significant character of the MRD is that drought and salinity intrusion normally appear in the dry season while flooding mostly occurs in the wet season. Nevertheless, it has a naturally suitable condition for agricultural cultivation and aquaculture, with approximately 26,000 km 2 of the total area as agricultural production land, corresponding to 65% of the total area of the region. The main water sources of irrigation are from the channel network and rainfall. Annually, it has contributed up to 55% of total production of paddy and 70% of total production of aquaculture. Therefore, the MRD plays an important role in the agricultural sector and food security of Vietnam.

Data Preparation
Main data used to compute TVDI are derived from the MOD09GQ and MOD021KM data products. The USGS provide these data products free of charge. MOD09GQ provides MODIS band 1-2 daily surface reflectance at 250 m resolution [50,51], where band1 covers a spectral range of 0.62-0.67 µm and band2 of 0.84-0.87 µm. These reflective data are used to compute NDVI. MOD021KM provides MODIS band 31-32 daily thermal emission at 1 km resolution [52,53]. These emissive data are used to compute LST. The observed period is in the dry season, from January to April, between 2001 and 2015.
To cover the MRD, we need two scenes at the locations of 'h28v07' and 'h28v08'. The couples of the two scenes were captured nearly at the same time. The images are referenced to the sinusoidal datum and stored in the hdf format. In this study, the collected data include 315 image couples in case of less than 10% of cloud cover during the observed period. Accordingly, the two scenes of image are merged to cover the whole region. Then, bands 1, 2, 31, and 32 are combined to store into one file. Moreover, the image dataset is converted to the WGS84 geographic coordinate system from the Sinusoidal datum. Eventually, the images are clipped keeping interior to the boundary of the MRD.

Data Processing
In this section, we describe the processing of the MODIS data to estimate the spatiotemporal trend of TVDI changes in the dry season between 2001 and 2015. For each image acquisition time, NDVI and LST are determined. Then, TVDI images are calculated based on NDVI and LST. Subsequently, a mean TVDI image is obtained for each dry season. Finally, for each pixel a temporal trend of TVDI changes in time series is estimated, using a linear regression. The resultant image is expected to be representative for the spatial pattern and temporal distribution of drought changes in the MRD in the dry season during the observed period.

NDVI
Based on Red and NIR reflective bands, NDVI is determined as in Equation (1). Here, ρ Red and ρ NIR are representative for the MDO09QG band 1 and 2 surface reflectance values respectively. Then, the NDVI image is reconstructed at the 1 km spatial resolution by the bi-linear interpolation. Figure 2a shows NDVI derived from the data on 12 February, 2015. (1)

LST
LST (K) is computed following Equation (2), based on the algorithm developed by Price (1984), and confirmed by Vazquez et al. (1997) [54,55]. Here, T 31 , T 32 (K) are the brightness temperatures obtained from band 31, band 32, respectively, and, ε 31 , ε 32 are the surface emissivity coefficients in band 31, band 32, respectively. In addition, the surface emissivity is calculated from NDVI, as Equations (4) and (5), apply the algorithm developed by Cihlar et al. (1997) [56]. Then, LST ( • C) is determined following Equation (6). Additionally, LST anomalies out of a range from 15 to 45 ( • C) are removed, conforming to the tropical weather in the study area. Figure 2b shows LST ( • C) derived from the data on 12 February, 2015.
The brightness temperature T b detected by a thermal sensor is determined by Planck's Equation (7).
where, L λ (Wm −2 sr −1 µm −1 ) is the spectral radiation, h = 6.62  2002) proposed an index, called TVDI, to monitor drought levels based on LST and NDVI [19]. The levels of drought were related to evapotranspiration, soil moisture (or LST), and vegetation coverage. LST will be smallest at surfaces corresponding to maximum evaporation due to saturated water, while it will increase to maximum at minimum evaporative surfaces due to very dry surfaces (with or without vegetation). Subsequently, TVDI is determined by Equation (8).
Here, LST ( • C) is the land surface temperature at the observed pixel. LST max and LST min depend on the NDVI value at the observed pixel, as Equations (9) and (10), corresponding to the temperatures at the dry edge and the wet edge, respectively. The parameters (a, b) and (c, d) are constant for each image acquisition time.
For estimation of the parameters (a, b), and (c, d), the linear regression was applied [58]. Firstly, the NDVI image was divided into 20 classes from 0 to 1 values with an interval of 0.05. For each NDVI class, a dataset of the LST values at the corresponding pixels was collected. Then, the maximum values LST max were chosen and the corresponding NDVI values NDVI max were determined, for an example see Table 1. Similarly, the minimum values LST min and NDVI min were chosen. Subsequently, the fitting lines were estimated from the datasets of (LST max , NDVI max ) and (LST min , NDVI min ), using the linear regression, as presented in the red line and blue line, see Figure 3. They are representative for the dry edge and the wet edge at the image acquisition time. Finally, the parameters (a, b) and (c, d) were derived from these two lines, for example equaling to (-11.5399, 37.9475) and (−0.2788, 22.3359) on 12 February, 2015, respectively.  According to Equation (8), Figure 4 shows the TVDI image in the MRD on 12 February, 2015. Here, TVDI values are classified to five levels, consisting of ranges from 0 to 0.2, 0.2 to 0.4, 0.4 to 0.6, 0.6 to 0.8, and 0.8 to 1, which are representative of drought levels such as wet, normal, light, moderate, and high, respectively [28,29].

Mean TVDI
For each dry season, mean and standard deviation (StdDev) images of TVDI are computed. For each pixel, the mean TVDI equals an average of at least six TVDI values, corresponding to at least six image acquisition times per season. Subsequently, the standard deviation TVDI is computed.
For example, Table 2 shows statistics of TVDI variations at the pixel of (10 • 06'N, 105 • 59'E), located at the Vinh Long Province. It consists of a number of TVDI images, mean and standard deviation of TVDI values for each dry season. As a result, this location shows a variation at moderate and high levels of drought in the dry season during the observed period.  (11), by the linear regression [58]. The trend is estimated if there are at least six mean TVDI values or seasons available between 2001 and 2015. The rate v of the trend is obtained by solving Equation (12) and the root mean square error (RMSE), as standard deviation of residuals, is also computed, as Equations (13) and (14). where, : the vector of the mean TVDI values per season.
x = x 0 v : the vector of parameters of the linear trend, offset x 0 and rate v.
: the design matrix, in which t i denotes the ith season.
RMSE: the root mean square error (RMSE), as standard deviation of residuals. e: the least-square residual vector. Figure 5 shows the temporal distribution of the mean TVDI values at the sampled pixel during the observed period. The result indicates that the location has an increase of drought change at a rate of +0.013 per year with RMSE of 0.05. Accordingly, a spatiotemporal pattern of drought changes in the dry season between 2001 and 2015 is estimated for the entire MRD.

The Mean TVDI Images
For each dry season, the mean TVDI image is computed and classified into the five levels of drought, and the total area of each drought level is determined in the entire MRD. Figure 6 presents the mean TVDI images in the dry season from 2001 to 2015. Most of areas present light and moderate levels of drought, corresponding to yellow and orange colors. A few of red small areas show a high level of drought appearing in urban and industrial areas. Inversely, light blue and blue areas appear in forest or wetland areas.
In details, Table 3 and Figure 7 present percentages of areas of drought levels in the MRD in the dry seasons during the observed period. The total area of the wet level occupies nearly zero percentage while the total area of the high level of drought varies from 1% to 2%. Obviously, the total area of the light and moderate levels mostly occurs in a range from 96% to 98%, occupying approximately 62% and 35%, respectively. Here, the area is calculated by the sum of the total pixels of the 1 × 1 km size.

The Spatiotemporal Pattern of Drought Changes
The spatiotemporal pattern of drought changes in the MRD in the dry season between 2001 and 2015 is shown in Figure 8a, in which each pixel presents a temporal trend of the mean TVDI values. Each pixel is classified into one of seven colored groups, based on its rate v of drought change, corresponding to the slope of the trend. Table 4 presents the total area of the seven v groups of drought change in the entire MRD. Accordingly, the total area of the pixels in green occupies 75% of the study area, meaning unchanged levels of drought. The pixels in blue-tone, occupying the total area of about 12.5%, present a negative trend of drought change. These areas can be wetter, occurring in forest, agricultural cultivation, and aquaculture. Most of them are in the border of the MRD, adjacent to the sea, including provinces Kien Giang, Tra Vinh, Soc Trang, and Bac Lieu. On the other hand, the red-tone pixels, also occupying the total area of about 12.5%, show a positive trend of drought change. Most of the areas with an increased trend appear in the center of the RMD, including provinces Dong Thap, Vinh Long, Can Tho, and Hau Giang. This means that these areas are prone to increased drought.    Table 4 indicates that sub-areas in Dong Thap, Vinh Long, Can Tho, and Hau Giang have an increased trend of drought change more the rest, occupying an area percentage of 30%, 45%, 29%, and 33%, respectively. Inversely, the sub-areas having a decreased trend are those in Bac Lieu, Kien Giang, Soc Trang, and Tra Vinh, with an area percentage of 27%, 28%, 23%, and 27%, respectively.

An Increase in Duration of the Dry Season
The climate of the MRD is tropical and dictated by two distinct seasons: the dry and wet seasons. The cool, dry air masses from north-eastern directions make the MRD much drier from January to March, which causes short droughts occurring sparsely over the entire region for two or three weeks. Meanwhile, due to the southwest Monsoon, the MRD mostly has more frequent and heavy rains from May to October, and then flooding usually occurs for a few weeks [59]. The average annual rainfall ranges from approximately 1000 to 2400 mm, and regularly decreases from southeast to northwest, e.g., approximately 2250 mm, 1600 mm, and 1490 mm at the meteorological stations in Ca Mau, Can Tho, and Dong Thap, respectively [43]. Due to reduction in rainfall and increase in hot weather duration during the dry seasons, the short droughts have become longer, denser, and more serious. Annually, a zero or very low rainfall period is normally experienced between January and March, but the average monthly rainfall of approximately 13 mm occurs sparsely over the MRD, as derived from the rainfall data at the meteorological stations in Can Tho and Dong Thap (data source: Vietnam Center of Hydro-Meteorological Data), see Figure 9. The increase in dry season duration is mostly caused by a rainfall shortage occurring in last December or April or both. For example, recently severe droughts in the MRD from 2002 to 2004 and 2015 to2016 were reported [44,59]. Based on the two monthly rainfall datasets at the Can Tho and the Dong Thap, the three-month SPI values through the end of March from 1985 to 2018 were determined, by applying the Standardized Precipitation Index Tool [60]. As a result, the SPI-3 values vary in a range from approximately −1.5 to 2.0, as presented in Figure 10. Accordingly, after the 1987 serious short drought, the MRD mostly experienced a normal state in the dry season with the SPI-3 range of (−1.0, +1.0), meaning no anomaly of drought or normal rainfall during the 10-year period. Nevertheless, positive and negative anomalies of drought have appeared extremely and alternatively since 2000, meaning that the precipitation regime is very dry or very wet, corresponding to the SPI-3 values of less than −1.5 or more than 1.5, respectively. For instance, short droughts occurred in 2002 and 2015 over the widespread MRD, corresponding to the SPI-3 of −1.5 at the Can Tho and −1.0 at the Dong Thap. Table 5 shows a comparison between the TVDI and SPI-3 indicators at the meteorological stations in Can Tho and Dong Thap during the observed period. The SPI-3 values present a large variation of drought anomalies from negative to positive while the TVDI values mostly illustrate moderate levels of drought, also see Figure 11. Typically, the SPI-3 indicator through the end of March only depends on accumulated rainfall in the dry season, from January to March, and the long-term rainfall from 1985 to 2018. Meanwhile, the TVDI indicator is affected by land cover characteristics and surface temperature, or evapotranspiration, during the dry season. Subsequently, the TVDI and SPI-3 values strongly indicate moderate and severe droughts when rainfall shortages appear, e.g., from 2002 to 2003, and 2014 to 2015.

Year
Can Tho Dong Thap

Spatiotemporal Pattern of LST Changes
Similarly, based on the dataset of the LST images, a spatiotemporal pattern of LST changes in the MRD in the dry seasons is estimated, by the linear regression, as shown in Figure 12. It means that for each pixel, a temporal trend of LST change is determined and its rate of LST change is described as the slope of a fitting line estimated from LST values in time-series. Each pixel is classified into one of five colored groups, based on its rate of the LST changes. As shown in Table 6, the total area of the pixels in yellow occupies about 50% of the study area, meaning that LST are mostly unchanged. The pixels in blue-tone, occupying the total area of about 5%, present a negative trend of the LST changes, while the pixels in red-tone, also occupying the total area of about 45%, show a positive trend of the LST changes. In general, LST in the MRD increases at an average rate of +0.1 • C per year in the dry season between 2000 and 2015. Subsequently, the pattern of LST changes in the MRD is positively correlated to the spatiotemporal pattern of drought changes during the observed period.

Drought Level Change vs. Land Use Change
Based on land use maps in the Southern Vietnam, investigated and created by the provincial governments every five years, Table 7 presents a summary of main land use areas belonging to the MRD in 2000 and 2015 [61]. Subsequently, the provinces, including Dong Thap, Vinh Long, Can Tho, and Hau Giang, have an obvious change in land use, as forest and land for cultivation of perennial trees decreased by 16.9, 83.6, and 86.0 km 2 while built-up land for residence and public transportation increased by 52.6, 49.3, and 82.0 km 2 , respectively. Inversely, there were also obvious land use change in the maritime provinces, e.g., Tra Vinh, Bac Lieu, Ca Mau, and Kien Giang, in the following three ways. Firstly, paddy land and land for cultivation of annual crops were converted to land for cultivation of perennial trees or aquaculture land. Secondly, forest land and/or land for cultivation of perennial trees were converted to aquaculture land nearby the coastal areas. Thirdly, mangrove forests were planted along most of the coastal areas. Changes in forest land are obviously recognized in the forest maps in the MRD in 2005 and 2015, as shown in Figure 13 [62].  Additionally, Pham et al. (2016) confirmed that agricultural land decreased, including land for cultivation of paddy and annual crops, perennial trees, and forest, while non-agricultural land increased, including residential land, land for construction, and land for public purposes, between 2005 and 2014 [62]. As a result, the two land use maps in the MRD were derived from the Landsat-7 scene on 20 January, 2005 and the Landsat-8 scene on 21 January, 2014, as presented in Figure 14. Then, the total areas of land use types were computed, as shown in Table 8. Accordingly, land for cultivation of paddy and annual crop, colored in yellow, occupying 54% of the total area in 2014, decreased by~1000 km 2 . Land for cultivation for perennial trees, most of fruits, colored in light orange, occupying 18.8%, decreased by~1400 km 2 . Forest land, colored in green, only occupying 3.1%, decreased by~625 km 2 , mostly occured in Long An, and Dong Thap. Inversely, the others, colored in gray, occupying 12.6%, including urban and rural residential land, land for construction of offices, industrial parks, and public transports, strongly increased by~3000 km 2 , clearly appeared in Long An, Tien Giang, Dong Thap, An Giang, Vinh Long, Can Tho, and Hau Giang. In summary, the sub-areas in the central region have a positive trend of drought change due to an increase of imperious surface, i.e., land for cultivation of perennial trees and forest have been converted to built-up land for residence and public transportation. Meanwhile, the sub-areas along the coastal region have a negative trend of drought change due to an increase of water and vegetated surface, e.g., mangrove and aquaculture land.

Conclusions
We have presented a state-of-the-art analysis of drought occurring in the MRD in the dry season from 2001 to 2015. In this case, the TVDI indicator is exploited to assess drought levels, using the daily MODIS data products. The spatiotemporal pattern of drought changes in the MRD during the observed period is estimated based on the mean TVDI images, using the linear regression. The result indicates that approximately 62% of the study area is at light level of drought while 34% are at moderate level in the dry season. Subsequently, approximately 12.5% of the total area have an increased trend of drought change, indicative of becoming drier, mostly occurring in the central sub-regions. Inversely, also approximately 12.5% of the total area are shown by negative trends of drought change, becoming wetter, appearing in several sub-areas in the maritime provinces. The results are expected to provide important information for the sustainable water resource management and planning in the MRD.