Geospatial Assessment of Agriculture Land Use Change in Relation to Paddy Production along the Punakha-Wangdue Valley over the Last Decade

1. Introduction

For agricultural sustainability and food security, the intricate relationship between land use dynamics and food production is critical. Paddy, a staple food for many countries worldwide is a lifeline for billions and is of huge economic importance (Muthayya et al., 2014). Land use change refers to the process of change in land cover associated with the expansion or contraction of the area of land used for a specific purpose such as pasture, cropland or urban (Davis et al., 2019). More than half of the global land surface has experienced at least some form of change throughout the last millennia (Luyssaert et al., 2014).

Land-use land-cover [LULC] changes have long caused and are attributed to global biodiversity loss, land degradation and climate change (Oliver & Morecroft, 2014). Ever since the unprecedented industrialisation during the nineteenth century, LULC changes occurred more frequently and to a greater extent through rapid deforestation, conversion of open mostly agricultural land into settlement in some areas and through broader exploitation of natural resources and intensification of agriculture (Foley et al., 2005). A great deal of contemporary urbanisation has been characterised as urban sprawl, which is a highly extensive form of land take for urban uses having environmentally detrimental effects (Nuissl & Siedentop, 2021).

Bhutan is no exception to increase in built-up areas due to urbanization which leads to encroachment towards areas under vegetation and agriculture (Yangchen et al., 2015).. According to Gyeltshen et al. (2022), the rice self-sufficiency ratio of Bhutan stands at 34.7-51.03% thereby relying the rest on import.

Punakha and Wangdue Phodrang are the major paddy producing districts in Bhutan with a combined share of 17% of the total paddy production of the country in 2021 (NSB, 2021). As per the Ministry of Works and Human Settlement [MoWHS] (2008), Wangdue Phodrang followed by Punakha has the highest developmental potential from the Central-Western regions of Bhutan. Such a potential poses a threat to the continuity of rice production in these districts. Therefore, it is crucial to understand the spatial trend of LULC changes and its impact on the paddy production scenario so that the state and potential risks of particularly encroachment of already scarce agricultural land can be highlighted and discussed.

2. Materials and Methods

2.1. Study Area

As shown in Figure 1, Punakha and Wangdue Phodrang districts are situated between 800 - 2500 metres above sea level in western-central Bhutan. The Punakha-Wangdue valley commonly experiences wet summers and dry winters due to its subtropical climate (Chhogyel & Bajgai, 2016). The favourable climatic and soil conditions of the districts allow for the growth of various cereals and fruits of which rice, wheat, mustard, potato and citrus are major cash crops cultivated (Punakha Dzongkhag Administration, 2019; Wangdue Phodrang Dzongkhag Administration, 2018).

2.2. Satellite Data Acquisition

Atmospheric correction of satellite data is an important step that needs to be conducted before it is used to reduce error caused by particles present in the atmosphere (Makarau et al., 2014). So, multi-spectral atmospherically corrected surface reflectance [SR] data of Landsat 7 Enhanced Thematic Mapper [ETM+] which were processed using Landsat Ecosystem Disturbance Adaptive Processing System [LEDAPS] algorithm version 3.4.0 for year 2005 and 2010 and atmospherically corrected surface reflectance data of Landsat 8 Operational Land Imager/Thermal Infrared Sensor [OLI/TIRS] using Land Surface Reflectance Code [LaSRC] for the year 2015 and 2020 was acquired respectively from the Google Earth Engine [GEE] data catalogue.

Figure 1. Study area.

Figure 1. Study area.

2.3. Cloud Masking

As cloud cover is an unavoidable problem which blocks the optical data of satellite images (Foga et al., 2017) cloud masking was also conducted through the bitmask function using the ‘QA_PIXEL’ band of the respective satellite data in GEE.

2.4. Image Classification

ArcGIS 10.8 was used to perform the LULC classification processes of the acquired images. A supervised classification method where pixels of images are compared and matched with known pixels of interest to assemble them and generate different classes was used to classify the images (Palaniswami et al., 2006). Fifty training sites were defined from the satellite image and merged together to create individual classes and then create a signature file to perform maximum likelihood classification. Five different classes were used to generate the LULC map of respective images namely Vegetation, Agricultural land, Snow cover, Built up and Barren area.

2.5. Accuracy Assessment of Classification

Accuracy assessment is a critical step in image classification in remote sensing to quantitively evaluate how the generated classes are assigned with correct real world land classes (Rwanga & Ndambuki, 2017). The accuracy of the image classification was determined through the kappa statistic as shown in equation 1 by generating 200 random points and comparing it with the basemap of ArcGIS (Rwanga & Ndambuki, 2017).

$$Kappacoefficient={\displaystyle \frac{\sum ki=1nii-\sum ki=1nii\left(GiCi\right)}{n2-\sum ki=1nii\left(GiCi\right)}}$$

where, i is the class number, n is the total number of classified pixels that are being compared to actual data, nii is the number of pixels belonging to the actual data class i, that were classified with a class i, Ci is the total number of classified pixels belonging to class i and Gi is the total number of actual data pixels belonging to class i.

3. Data Analysis

3.1. Change Detection of LULC

The area statistics of all the classes of classified images were calculated using the raster calculator of respective years in ArcGIS. The rate of change of land usage from 2005 to 2020 use was determined by categorizing the study years into three time periods 2005-2010, 2010-2015 and 2015-2020 as a percentage of decrease or increase from the previous year as shown in equation 2 (Kafi et al., 2014).

$$Changeinarea\left(\%\right)={\displaystyle \frac{TotalareainYii-TotalAreaYi}{TotalareaYii}}\times 100\dots \dots \mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}2$$

where Yi is the first year and Yii is the second year.

3.2. Relative Conversion Index [RCI] of Agricultural Land

The intensity of agriculture land conversion in Punakha-Wangdue valley to other land uses was quantified using the Relative Conversion Index (equation 3) as described by Susilo (2016).

$$RCI={\displaystyle \frac{RDk}{RTk}}\dots \dots \dots \dots \mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}3$$

$$RDk={\displaystyle \frac{ADk}{AD}}\dots \dots \dots \dots (3.1)$$

$$RTk={\displaystyle \frac{ATk}{AT}}\dots \dots \dots \dots (3.2)$$

where RDk is the proportion of conversion of agricultural land in each sub-district (Gewog) (equation 3.1) and RTk is the proportion of agricultural land conversion in the entire study area (equation 3.2). The intensity of agricultural land conversion was also visually represented on a map by categorizing the RCI values as Low, Moderate and High (Figure 4).

3.3. Change in Paddy Production

Data of total paddy production from 2005 to 2021 for Punakha and Wangdue Phodrang respectively was extracted from annual agriculture statistics publications and was used to calculate the change in paddy production as a percentage from a previous year to a given year.

$$Changeinproduction\left(\mathrm{\%}\right)={\displaystyle \frac{TotalproductioninYii-TotalproductionYi}{TotalproductionYi}}\times 100\dots .\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}4$$

where Yi is the first year and Yii is the second year.

3.4. Calculation of Rice Self-Sufficiency Ratio

The rice self-sufficiency ratio for Punakha and Wangdue Phodrang was carried out using equation 5 as described by Gyeltshen et al. (2022) for 2016 to 2021.

$$Self-sufficiencyratio={\displaystyle \frac{DRP}{TRC}}\times 100\dots \dots \mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}5$$

where, DRP is the domestic rice production and TRC is the total rice consumption.

The domestic rice production was calculated by obtaining total paddy production data for Punakha and Wangdue Phodrang respectively from 2016 to 2021 from NSB (2021) and adjusted for rice milling recovery of 65% (Chhogyel et al., 2015; Gyeltshen et al., 2022).

The total rice consumption was determined using equation 6 as described by Firmansyah et al. (2021) for 2016 to 2021.

$$Totalriceconsumption=percapitaconsumption\times population\dots \dots \mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}6$$

Since there was no official rice per capita consumption figure for Bhutan, the average national rice consumption per capita was taken as 147.7 kg as calculated by Gyeltshen et al. (2022). The regional population data for Punakha and Wangdue districts from 2016 to 2021 were obtained from NSB (2019).

Figure 2. Gain/Loss percent of individual classes between 2005-2020.

Figure 2. Gain/Loss percent of individual classes between 2005-2020.

4. Results and Discussion

4.1. Accuracy Assessment

The kappa coefficients obtained for respective years are shown in Table 1.According to McHug (2012) kappa coefficient between 0.60 to 0.79 is substantial and the kappa coefficient between 0.80 to 0.90 represents a strong level of agreement. Since the kappa coefficient obtained for respective LULC maps falls between 0.60 and 0.90, the classification accuracy can be considered to have a substantial to strong level of agreement.

4.2. Change in Land Use Land Cover Analysis

There was an overall increase in vegetation (Figure 2) where similar study conducted by Gilani et al. (2015) also found that there was a net increase in vegetation cover in Bhutan with a gain in forest cover far outweighing the loss during the entire twenty-year study period from 1990 to 2010 in which major gains were observed in Bumthang and Wangdue Phodrang districts. The global forest resource assessment 2020 report Bhutan published by FAO (2020) also reported an overall increasing trend of vegetation in the country from 1990 to 2020 which aligned with Bhutan’s reforestation effort and activities carried out during the period. However, a slight decrease in vegetation in recent years from 2015 to 2020 in the region could be attributed to anthropogenic factors such as an increase in built-up areas and other developmental activities causing deforestation and thus reduction in vegetation (Dorji et al., 2022). Multiple forest fires that occurred in 2018 in Wangdue Phodrang could also have contributed to the loss of vegetation (Zangmo, 2018). Regardless, changes in forest cover remain fairly stable due to the country's conservation efforts and reforestation schemes as directed by its development philosophy (Bruggeman et al., 2016).

Agricultural land on the other hand reduced by 5.1% (Figure 2). The findings are in line with Yangchen et al. (2015) where the authors found a significant reduction in agriculture area in the country from 14% to 8% from 2000 to 2013. Another study conducted by Nepal (2024) to see the LULC change in Bhutan from 1930-2020 also points to the reduction in agricultural area citing the influence of urbanization and shift in land-use practices. The development of Bajothang town in 2008 required the conversion of 15 ha of fairly flat agricultural wetland which added to the declining trend of agriculture area in the region (Wangzom & Zangmo, 2018). According to Tenzin et al. (2019), challenges posed by climate change such as erratic rainfall patterns, windstorms and reduced water availability are a threat, particularly for Bhutanese agriculture systems.

Bhutan has seen a gradual increase in built-up areas over the decades as major planning for infrastructure development started in the 1960s (Wangchuk & Turner, 2019). The LULC analysis reveals that there was a rise in builit-up area by 4.8% in the region (Fig 1). As growing population leads to more and rapid development pressure for basic amenities, transportation facilities, entertainment and others infrastructure (Tobgye & Piyathamrongchai, 2019) the growth of urban population from 2,292 to 7,046 and 7,522 to 9,880 in Punakha and Wangdue Phodrang respectively could have driven the incremental growth of built up in the region. The reallocation and development of the new Bajothang town phase one also could have contributed to the rise of built up areas in the region (Wangzom & Zangmo, 2018). Construction of major hydropower projects such as Punatshangchu Hydroelectric Project I [PHP-I] that started in 2008 and PHP-II in 2010 in the region have also required the development of offices and labour settlements causing the built up areas to expand (Tshering, 2012).

Barren area in the region also increased (Figure 2) over the study period which according to Wang et al. (2018) could be attributed to the increased fallowing of lands that are adjacent to urban areas in anticipation of future town growth. Exposure of bare soils due to excavation and debris dumping activities carried out during huge construction projects such as PHP-I and PHP-II could have added to the barren area in the later year of the study period in the region (Wangdi, 2014).

The findings of the overall declining trend of snow cover in the region are in accordance with the findings of Gurung et al. (2011) where average Snow Cover Area [SCA] for the period between 2002 to 2010 was found to be decreasing in the country. Similar trends were also observed by Mahagaonkar et al. (2017) where they suggest evidence of extensive glacial retreats due to warming temperatures in the Himalayas of Bhutan. According to Slatyer et al. (2022), the rise in global temperature is one of the major factors responsible for widespread decline in snow cover regimes around the world. So, the decline in snow cover in the region could be attributed to rising air temperature over the study period in the region (NCHM, 2019).

Figure 3. LULC transition map between 2005-2010, 2010-2015, 2015-2020.

Figure 3. LULC transition map between 2005-2010, 2010-2015, 2015-2020.

Table 2. Change in area of individual classes between 2005-2010 and 2010-2015.

Class	2005-2010			2010-2015
	Area (ha)		Change (%)	Area (ha)		Change (%)
	From	To		From	To
Vegetation	348983.3	355595.1	1.9%	355595.1	355906.6	0.1%
Agriculture	12660.6	12539.5	-1.0%	12539.5	12281.7	-2.1%
Water Body	2114.2	2218.0	4.9%	2218.0	2291.9	3.3%
Snow Cover	30728.1	27235.3	-11.4%	27235.3	28210.4	3.6%
Built Up	1155.5	1169.0	1.2%	1169.0	1178.9	0.8%
Barren Area	119652.5	116537.2	-2.6%	116537.2	115424.6	-1.0%

Table 2. Change in area of individual classes between 2005-2010 and 2010-2015.

Class	2005-2010			2010-2015
	Area (ha)		Change (%)	Area (ha)		Change (%)
	From	To		From	To
Vegetation	348983.3	355595.1	1.9%	355595.1	355906.6	0.1%
Agriculture	12660.6	12539.5	-1.0%	12539.5	12281.7	-2.1%
Water Body	2114.2	2218.0	4.9%	2218.0	2291.9	3.3%
Snow Cover	30728.1	27235.3	-11.4%	27235.3	28210.4	3.6%
Built Up	1155.5	1169.0	1.2%	1169.0	1178.9	0.8%
Barren Area	119652.5	116537.2	-2.6%	116537.2	115424.6	-1.0%

Table 3. Change in area of individual classes between 2015-2020 and 2005-2020.

Class	2015-2020			2005-2020
	Area (ha)		Change (%)	Area (ha)		Change (%)
	From	To		From	To
Vegetation	355906.6	352366.0	-1.0%	348983.3	352366.0	1.0%
Agriculture	12281.7	12023.5	-2.1%	12660.6	12023.5	-5.0%
Water Body	2291.9	2208.1	-3.7%	2114.2	2208.1	4.4%
Snow Cover	28210.4	27445.7	-2.7%	30728.1	27445.7	-10.7%
Built Up	1178.9	1211.7	2.8%	1155.5	1211.7	4.9%
Barren Area	115424.6	120078.4	4.0%	119652.5	120078.4	0.4%

Table 3. Change in area of individual classes between 2015-2020 and 2005-2020.

Class	2015-2020			2005-2020
	Area (ha)		Change (%)	Area (ha)		Change (%)
	From	To		From	To
Vegetation	355906.6	352366.0	-1.0%	348983.3	352366.0	1.0%
Agriculture	12281.7	12023.5	-2.1%	12660.6	12023.5	-5.0%
Water Body	2291.9	2208.1	-3.7%	2114.2	2208.1	4.4%
Snow Cover	28210.4	27445.7	-2.7%	30728.1	27445.7	-10.7%
Built Up	1178.9	1211.7	2.8%	1155.5	1211.7	4.9%
Barren Area	115424.6	120078.4	4.0%	119652.5	120078.4	0.4%

4.3. Relative Conversion Index Analysis of Agriculture Land

The areas with higher intensity of agricultural land conversion are spatially situated in the central regions of the study area whereas areas with lower intensity of agricultural land conversion are situated at the northern and eastern outskirts of the Punakha-Wangdue region (Figure 4). According to Chen et al. (2021), urbanisation influences the information entropy of land use structure disturbing the orderly natural landscape causing the land use intensity and transition to rise as the area transitions into a semi-artificial landscape before becoming fully artificial.

Figure 4. Intensity of agriculture land conversion in Punakha and Wangdue Phodrang from 2005 to 2020.

Figure 4. Intensity of agriculture land conversion in Punakha and Wangdue Phodrang from 2005 to 2020.

This phenomenon explains the concentration of Gewogs having high agricultural land conversion intensities at the central lower region as there are numerous clusters of settlements present at the location with key towns that are Bajothang and Khuruthang at Thoedtsho and Guma Gewogs respectively. Similarly, Bapisa Gewog encompassing key roadside stop point town Lobesa which connects people travelling via road from northern, southern and eastern Bhutan has also experienced significant conversions. Such influences were also observed by Naab et al. (2013) where agricultural lands were repurposed in peri-urban areas.

4.4. Change in Paddy Production

Paddy production in the region shows a fluctuating trend over the years from 2005 to 2021. There was a net reduction of paddy production in Wangdue Phodrang districts by 7.2% between 2005 to 2021. The district recorded the highest paddy production in 2014 with 9173 metric tonnes while the lowest was in 2011 with 2783 MT. Punakha district also had a decrease in overall paddy production by 5.9% compared to 2005 in 2021. In 2018, the district achieved its peak production of 16,389 MT of paddy, whereas the lowest production was observed in 2021, amounting to 6,509 MT (Figure 5).

Similar trend of initial increase in paddy production and decrease in recent years was also observed Gyeltshen et al. (2022). The initial positive trend of paddy production could be due to efforts taken by the government and agriculture research centres to release improved and high yielding rice varieties which in turn increased the productivity (per area production) (Tshering et al., 2021) in the districts as Chhogyel and Bajgai (2016) reported that the adoption of modern rice varieties stood at 68% and 62% in Wangdue Phodrang and Punakha respectively. However, reduced production in recent years could be attributed to decrease in harvest area due to increased fallowing (Ngawang, 2018), human wildlife conflict (Wangchuk et al., 2023), labour shortage and high production costs (Gyem et al., 2018; Tashi & Dendup, 2020), weeds and diseases (Tshewang et al., 2016), climate change (Chhogyel & Kumar, 2018) and inadequate water sources and irrigation infrastructure (Kusters & Wangdi, 2013).

4.5. Impact of LULC Change and Paddy Production on Rice Self-Sufficiency

From the rice self-sufficiency estimates drawn by calculating total rice consumption and production within the districts (Punakha and Wangdue Phodrang) with adjustments taken for rice milling recovery of 65%, it showed a declining trend for both Punakha (Figure 6) and Wangdue Phodrang (Figure 6) from 2016 to 2021. Over this period, fluctuations in production and consumption dynamics are evident. In 2016 and 2017, Punakha recorded significant surpluses, with production surpassing consumption by 127% and 120%, respectively, indicative of robust agricultural output exceeding local demands. However, subsequent years saw varying trends, with instances of both surplus and deficiency. Notably, in 2021, a deficiency of 6% was recorded, concerning the state of food security within the district (Figure 6).

Similarly, the results reveal a consistent pattern of deficiencies in rice production compared to consumption throughout 2016 to 2021 in Wangdue Phodrang (Figure 6). In 2016, the deficit stood at 17%, indicating a shortfall in meeting consumption demands. This trend persisted in subsequent years, with deficits of 8%, 32%, 42%, 40%, and 57% recorded in 2017, 2018, 2019, 2020, and 2021, respectively. These figures signify a widening gap between rice production and consumption, culminating in a significant deficiency by the year 2021. Such persistent deficiencies raise substantial concerns regarding food security within Wangdue Phodrang, necessitating urgent attention and strategic interventions to address rice production.

Bhutan has been striving towards attaining rice self-sufficiency though figures show otherwise with increasing rice imports with rice self-sufficiency ratios [SSR] calculated well below 50% between 2005 and 2019 (Gyeltshen et al., 2022). According to Biswajit et al. (2013), several other neighbouring South Asian countries such as Nepal, Bangladesh and Sri Lanka are also not rice self-sufficient except for India and Pakistan. Similarly, most African countries are also reliant on imports to meet local rice demand (Van Oort et al., 2015).

A variety of factors are responsible for influencing the paddy production and ultimately the rice SSR of a country. In a study conducted by Dawe (2013) to assess geographic determinants of rice self-sufficiency, he found that available land suitable for rice cultivation those are flat with enough water resources influences per capita production and harvest area thus increasing the rice self-sufficiency compared to those with challenging topography. This highlights the impacts that land use changes have on paddy production and ultimately rice self-sufficiency for particularly those countries with difficult terrain where suitable land is already limited as in Bhutan and Nepal (Bhandari et al., 2023). Other factors include nutrient deficiencies, poor harvest, rising cost of production, depleting water resources and environmental causes (Kumar et al., 2021). Challenges such as lack of transportation facilities, inefficient irrigation systems, and inadequate land for cultivation are also listed as core constraints for paddy cultivation in rural areas of Malaysia (Echoh et al., 2017).

In the Bhutanese context, pertinent rice production challenges apart from limited arable agricultural area are land use change where social and economic development is causing the loss of prime wetland and the preference of farmers to grow more economically competitive and less resource intensive crops instead (NSB, 2021), rural-urban migration causing labour shortages driving costs (Yoezer, 2019), vulnerability of Bhutanese rice farming systems to climate change impacts as a major portion of rice growers are dependent on monsoon and spring water for irrigation (Chhogyel & Kumar, 2018) and price volatility due to limited domestic and international trade options (Chhogyel et al., 2015; Gyem et al., 2018).

All these findings and evidences solidify the direct relation of land use change on paddy production within the districts as the region experienced declining production with increasing land use change be it fallowing, conversion to other uses or shift in agriculture systems where alternative crops are cultivated instead of paddy over the years. One key take from this study is the central role land use change have on paddy production as it encompasses all other influencing factors thus determining the paddy production (Figure 7) as the harvest area is strongly correlated to production (Yasar & Siwar, 2016). With the districts and country as a whole experiencing change in land use of which fallowing and conversion are dominant (Ngawang, 2018) it is a growing threat to agricultural production including paddy with findings from this study pointing its influence on production and its impact on rice self-sufficiency over the years.

5. Conclusions

The geospatial assessment revealed that major land use land cover changes happened in Punakha and Wangdue Phodrang districts from 2005 to 2020. Relative conversion index analysis of agriculture land indicated that higher conversion intensity was spatially located in the central areas of the region while lower intensity was observed at the outskirt areas of the region. Yearly estimates of the rice self-sufficiency ratio within the districts reveal an increasing deficit trend with a rice surplus of 127% in 2016 falling to a deficit of 6% in 2021 for Punakha while previously rice deficit of 17% in 2016 declined further to 57% in 2021 for Wangdue. Such trends though influenced by the interaction of various biophysical, social, economic and political factors, demand urgent intervention to reduce impacts.