Environmental drivers of spatio-temporal dynamics in floodplain vegetation: grasslands as habitat for megafauna in Bardia National Park (Nepal)

Disturbance-dependent grasslands, often associated with hydromorphological and fire dynamics, are threatened, especially in subtropical climates. In the Nepalese and Indian Terai Arc Landscape at the foot of the Himalayas, natural and cultural grasslands serve a viable role for greater one-horned rhinoceros (Rhinoceros unicornis) and for grazers that form prey of the Royal Bengal tiger (Panthera tigris). The grasslands are vulnerable to encroachment of forest. We aimed to establish the effects of environmental drivers, in particular river discharge, river channel dynamics, precipitation and forest fires, on the spatio-temporal dynamics of these grasslands. The study area is the floodplain of the eastern branch of the Karnali River and adjacent western part of Bardia National Park. We created annual time series (1993–2019) of land cover with the use of field data, remotely sensed LANDSAT imagery and a supervised classification model. Additionally, we analysed the pattern of grassland patches and aerial photographs of 1964. Between 1964 and 2019, grassland patches decreased in abundance and size due to encroachment of forest. Outside the floodplain, conversion of grassland to bare substrate coincides with extreme precipitation events. Within the floodplain, conversion of grassland to bare substrate correlates with the magnitude of the annual peak discharge of the bifurcated Karnali River. Since 2009, however, this correlation is absent due to a shift of the main discharge channel to the western branch of the Karnali River. Consequently, alluvial tall grasslands (Saccharum spontaneum dominant) have vastly expanded between 2009 and 2019. Because the hydromorphological processes in the floodplain have become more static, other sources of disturbances – local flooding of ephemeral streams, anthropogenic maintenance, grazing and fires – are more paramount to prevent encroachment of grasslands. Altogether, our findings underscore that a change in the environmental drivers impact the surface area and heterogeneity of grassland patches in the landscape, which can lead to cascading effects for the grassland-dependent megafauna.

1 Introduction

Grasslands occur in diverse environmental conditions that can be dry, wet, hot, cold, productive, barren, dynamic and static (Brown and Makings, 2014). Their global areal extent has decreased by 40 % since the Industrial Era (White et al., 2000) due to cultivation, urbanization and increased forest growth, causing loss of biodiversity (White et al., 2000; Veldman et al., 2015). Grasslands are home to most of the extant mammalian megafauna on Earth and provide food for herbivores, which in turn sustain their predators. This is for example the case in the African (sub)tropical grasslands and in the mountainous and colder grassland habitats in Mongolia and Nepal. Here, grassland-dependent blue sheep (Pseudois nayaur), Siberian ibex (Capra sibirica) and Tibetan argali (Ovis ammon) serve as principle prey species of the snow leopard (Panthera unica) and Himalayan wolf (Canis lupus chanco) (Oli, 1996; Lyngdoh et al., 2014; Chetri et al., 2017; Rovero et al., 2020; Filla et al., 2021). At the foot of mountain ranges, warmer and wetter conditions are present due to lower altitudes and orographic lift of seasonal moist winds, causing the constraints for forest growth to be lifted (Staver et al., 2011; Hirota et al., 2011). The relationship between prey, predator and grassland habitat is also present in these regions more suited for forest growth: grassland patches in forests at the foot of mountainous regions in South Asia sustain Royal Bengal tigers (Panthera tigris) and their prey. This is for example the case in the Western Ghats, southern India (Sankaran, 2009), in the Duars, eastern India and Bhutan, and in the Terai Arc Landscape (TAL), a subtropical belt at the foot of the Indian and Nepalese Himalayas (Harihar et al., 2014, Irengbam et al., 2017). The TAL has also seen a reduction of grassland in protected areas (DNPWC, 2016).

Grasslands in the TAL can broadly be subdivided into (1) natural alluvial tall grassland, located near active and in abandoned river channels, and (2) a mosaic of short grassland (phantas) and mixed tall grassland, located further away from the river channels. The alluvial tall grasslands serve as essential habitat for the greater one-horned rhino (Rhinoceros unicornis) and hog deer (Hyelaphus porcinus) (Jnawali and Wegge, 2000; Odden et al., 2005; Pradhan et al., 2008; Thapa et al., 2013; DNPWC and DFSC, 2018), and their numbers would be greatly affected if changes occur in composition and areal extent of these grasslands (Peet, 1997; Odden et al., 2005). Tiger density is highest near grasslands and the short grasslands in particular provide nutritious feeding grounds for chital deer (Axis axis), their most important prey (DNPWC and DFSC, 2018). In addition to the surface area, the spatial pattern of grassland in the landscape is of importance: habitat heterogeneity is positively correlated with the abundance and habitat preference of chital deer in Chitwan National Park (Nepal) (Bhattarai and Kindlmann, 2012).

A number of grasslands in the TAL are considered to be maintained via physical disturbances that prevent succession to mature vegetational stages (Seidensticker, 1976; Dinerstein, 1979a; Lehmkuhl, 1989, 1994; Peet et al., 1999b). In general, the pattern of vegetation is regulated by autogenic and allogenic factors (Tilman, 1988), of which the latter are either environmental or anthropogenic. Allogenically, the community dynamics can be controlled via disturbances, of which the intensity and spatial extent have been identified as principal factors (Turner et al., 1998). The natural disturbances recorded in the TAL include hydromorphological processes, fires, climatic variables (precipitation, temperature, droughts) and herbivory, whereas the anthropogenic disturbances include cutting of grasses and trees for maintenance and resource collection, clearance by fires, livestock grazing, and, in the past, cultivation and creation of settlements (Seidensticker, 1976; Dinerstein, 1979a; Lehmkuhl, 1989, 1994). Changes in disturbance regimes impact other disturbance regimes (Van Langevelde et al., 2003; Trauernicht et al., 2013; Orem and Pelletier, 2015), the vegetation pattern and thus the ecosystem as a whole (Parr et al., 2012). Conversion of disturbance-dependent grassland to forest can occur within several years when changes occur in rainfall, fire or herbivory (Bond, 2008). In Bardia National Park (BNP), a conservation area in the TAL located near the Karnali megafan in western Nepal, grasslands have been declining due to encroachment of forest (Peet et al., 1999a; Jnawali and Wegge, 2000; Odden, 2007). From a hydrological perspective, there are indications that a major change occurred during the monsoon in 2009 (Sinclair et al., 2017): near the apex of the highly dynamic braided river system of the Karnali River, the dominant discharge channel relocated from its eastern branch (the Geruwa River that borders BNP) to the western branch (Kauriala River). The susceptibility of land cover to changes in disturbance regimes not only holds for the grasslands in the TAL but also for other ecosystems where climate is conducive to forest growth. In the Brazilian Pantanal, where similar grasslands are found at the foot of the Brazilian Highlands, extreme hydrologic conditions have been explicitly identified with remote sensing as a driving factor of the pattern of vegetation communities (Arieira et al., 2011). The Taquari megafan in the Pantanal experienced a shift in hydrological conditions due to an avulsion, increasing susceptibility to deforestation and fire (Louzada et al., 2020).

Wide-ranging progress has been made in understanding the drivers of vegetation dynamics through space and time (Hupp and Osterkamp, 1996; Lorenz et al., 1997; Turner et al., 1998; Corenblit et al., 2007; Vesipa et al., 2017). The composition and spatial distribution of riparian vegetation communities are heavily influenced by flood frequency, duration, intensity (Hupp and Osterkamp, 1996; Poff et al., 1997) and timing (Newbold and Mountford, 1997). Consequences of high discharges for vegetation can be physical damage and uprooting, anoxia due to prolonged inundation and burial by fresh alluvium. Low discharges affect vegetation with drought stress when the ground water table becomes too deep. These processes are integrated in concepts such as “fluvial biogeomorphic succession”, describing the reciprocal interactions between fluvial landforms and vegetation (Gunderson, 2002; Corenblit et al., 2007), “flood pulse concept”, describing lateral connectivity and recognizing that flow and its variability are the main drivers of ecological processes in floodplains (FPs, Junk et al., 1989), and “shifting-mosaic steady state” for vegetation, which specifies the proportion of distinct successional land cover in floodplains (Arscott et al., 2002; Kollmann et al., 1999). Additionally, fire dynamics (Flannigan and Wotton, 2001; Hoetzel et al., 2013; Iglesias et al., 2014), herbivory (Owen-Smith, 1988; Van Langevelde et al., 2003; Allred et al., 2011) and anthropogenic interferences (Wang et al., 2015) have been studied. The interrelations of environmental drivers are complex (Lehmann et al., 2014), and there is not yet a complete understanding. It is, therefore, of interest to study vegetation and its drivers at a high temporal resolution over a long time span and in its historic context. Literature in this domain is particularly limited to a small number of geographic regions (Dufour et al., 2019), and the subtropical grasslands in conservation areas in the TAL such as BNP are underrepresented. It is especially vital to understand the functioning of such a valuable ecosystem in relation to possible changes in drivers (such as the indicated redistribution of discharge in 2009) in order to provide insight for nature management strategies.

A promising method for studying land cover dynamics and its drivers is the use of Earth observation techniques, which enable diachronic analysis (Lallias-Tacon et al., 2017; Dufour et al., 2019; Harezlak et al., 2020; Solins et al., 2017; Basumatary et al., 2021; Louzada et al., 2020; Van Iersel, 2020; Corenblit et al., 2010; Dufour et al., 2015). Vegetation in and near floodplains can be mapped at the following scales and characteristics: type (Alaibakhsh et al., 2017), species composition (Rapinel et al., 2019; Plakman et al., 2020), physiology (Wagner-Lücker et al., 2013) and structure (Straatsma and Baptist, 2008; Jalonen et al., 2015). Regionally, remote sensing has been successfully used to monitor similar grasslands (Sarma et al, 2008; Biswas, 2010; Thapa, 2011), fires (Takahata et al., 2010; Ghimire et al., 2014) and channel migration of the Karnali River (Rakhal et al., 2021).

We aim to quantify spatio-temporal land cover change and establish relationships with the environmental drivers with focus on the ecologically important grasslands. We address both grasslands within the floodplain of the Geruwa River (the eastern branch of the Karnali River) and the grassland patches further away. Specifically, we give answers to the following sub-questions: (a) what is the spatio-temporal pattern of land cover, (b) what is the spatio-temporal variation in environmental drivers during the last 3 decades, and (c) what are the effects of environmental drivers on land cover change and what are possible mechanisms that explain these effects? We do so by mapping the historic development including recent annual dynamics of land cover with the use of field and remotely sensed data (LANDSAT). We then relate the land cover dynamics to abrupt changes in hydrological, meteorological and forest fire variables. The hypothesis is that between 1964 and 2019 a transition occurred from grasslands to later successional stages. The discharge distribution during peak discharges at the upstream bifurcation of the Karnali River would be a relevant driver for the vegetation pattern and associated habitats of megafauna.

2 Study area

2.1 Bardia National Park and Karnali River

Bardia National Park is located in the southwest of Nepal (Fig. 1) and was established as conservation area on 8 March 1976. The study area (135 km² of the now in total 968 km² of BNP) is delimited by (1) the Siwalik Hills in the north, where ephemeral streams originate that are tributaries to the Karnali River, (2) a buffer zone and settlements in the southeast, and (3) the Geruwa Branch of the Karnali River in the west, which is a highly dynamic braided river system (Sinclair et al., 2017) located on the Karnali megafan (USAID, 2018). The Karnali River rises from Mt. Kailash on the Tibetan Plateau as one of the three biggest rivers draining Nepal. A few kilometres south of the Siwalik Hills near the town of Chisapani (Fig. 1), the Karnali River bifurcates into the Geruwa River (eastern branch) and Kauriala River (western branch), joining together in India as a tributary to the Ganges River. The subtropical climate in the region encompasses the monsoon, cool and dry post-monsoon, and a hot and dry pre-monsoon (USAID, 2018). Temperatures reach a maximum of 45 ^∘C in the hot season and fall to 10 ^∘C in January (Bolton, 1976). The mean annual precipitation is 1560 mm (DHM, 2017).

Figure 1 Delineation of Bardia National Park (dark green) and the study area (red) in and near the floodplain of the Karnali River. The megafan of the Karnali River lies between the two branches (the Kauriala and the Geruwa River) at the west side of Bardia National Park. This is a false colour composite (RGB = 656 of LANDSAT 8 imagery, 2020).

2.2 Drivers of grasslands in the Terai Arc Landscape and Bardia NP

Drivers in the form of physical disturbances vary in type, temporal components (frequency, duration, and timing), their spatial component (location, extent) and their intensity (White, 1985; Turner et al., 1998). In the TAL, both natural and anthropogenic disturbances are regarded as important drivers for grasslands in shaping the spatial and temporal patterns of vegetation. In BNP, fluvial processes are regarded as being important at the landscape level. Flooding of the Karnali River and associated clearance of vegetation, together with deposition of fresh alluvium, have been observed (Seidensticker, 1976; Dinerstein, 1979a; Lehmkuhl, 1989, 1994; Peet et al., 1999b). Fires occur in BNP on a yearly basis and can be natural as well as anthropogenic. Forest fires occur most frequently in the Sal forest (Shorea robusta dominated) and in the months April and May (Ghimire et al., 2014). The phantas are annually set on fire by park staff (Dinerstein, 1979b; Peet, 1997). Appendix A provides an overview from the literature of the activities by locals and park staff (see also Bhatta, 2000; Thing et al., 2017). After establishment as a conservation area in 1976, free access for locals, who extensively relied on BNP for resources, was refused. Later on, active habitat management was initiated on grasslands to compensate for the decreased disturbances. Permits that allow resource collection are needed for locals to enter the park for a certain number of days. Since 1983, the days of access and the number of permits issued increased until 1994 and 1999, respectively, and hereafter decreased drastically (Bolton, 1976; Bhatta, 2000). Whether the amount of permits in a given year and the days of access are proportional to the amount of removed biomass is not known (Bhatta, 2000). Before and around 1976, extensive grazing was recorded by cattle and wild ungulates in the southwest of the study area. The very short grazing lawns that are present in BNP are a result from patch-selective grazing by ungulates after cutting or burning (Thapa et al., 2021). Large herbivores can strongly influence the pattern and community structure of vegetation (Owen-Smith, 1988; Dinerstein, 1992). In the park, Asian elephant (Elephas maximus) (Ram and Acharya, 2020; Shrestha and Shrestha, 2021) and greater one-horned rhinoceros (rhino census in 2021) were virtually absent, but their numbers increased to 113 and 38 in 2021, respectively. Ungulates have been censused a number of times in BNP (Dinerstein, 1980; Wegge and Storaas, 2009; Dhakal et al., 2014; Karki et al., 2016; Kral et al., 2017; DNPWC and DFSC, 2018, 2022). Some uncertainty is present around the numbers of ungulates, which have been debated in the literature (van Lunenburg et al., 2017; Wegge et al., 2019). In BNP, the recorded densities of chital deer (the most abundant mid-sized grazer) are 53.99 (standard error, SE, 10.29) (Dhakal et al., 2014), 29.3 (SE 4.3) (Karki et al., 2016), 56.44 (SE 5.75) (DNPWC and DFSC, 2018) and 44.05 (SE 6.3) (DNPWC and DFSC, 2022).

Figure 2 Flowchart showing steps from data sources to analysis.

2.3 Vegetation types and successional stages

Dinerstein (1979a) described six major vegetation types in the park, later modified into seven types (Jnawali and Wegge, 1993; Pokheral, 1993): alluvial tall grasslands (Saccharum spontaneum dominant), mixed tall grasslands (wooded grasslands), short grasslands (<2 m, previously cultivated fields or “phantas”), khair–sissoo forest, moist riverine forest, mixed hardwood forest and Sal forest. More detailed grassland associations are clarified in Appendix B (Lehmkuhl, 1989, 1994; Peet et al., 1999a; Dinerstein, 1979a). The alluvial tall grasslands, which can be up to 4 m in height, consist of pioneer species that quickly cover bare alluvium after disturbance events. With no disturbances present, the line of succession is from alluvial tall grasslands via khair–sissoo forest to other riverine forest types, mixed hardwood forest or the climax Sal forest. Khair–sissoo forest is prevalent along the river courses and on floodplain islands in the Terai due to its ability to withstand floods. Sal forest consists of predominantly Shorea robusta trees that are fire resistant and is the most dominant tree species in BNP. Disturbances prevent succession of alluvial grasslands to forest, turning them into mixed tall grasslands. If disturbances continue, short (<2 m) and open grasslands may appear (locally called “phantas”), and their location is often more distal from stream channels (Dinerstein, 1979a; Lehmkuhl, 2000). A number of these grasslands are considered to owe their existence to human interferences, such as the Lamkauli, Baghaura and Khauraha phantas (Fig. 1) (Pokheral, 1993; Dinerstein, 1979b; Peet et al., 1999a). The largest phantas present are Lamkauli, Baghaura and Khauraha (Fig. 1). Nowadays, the phantas are predominantly vegetated with Imperator cylindrica, Vitiveria zizanioides and Desmostachyia bipinnata (Peet et al., 1999b). Within the mosaic of short and mixed tall grasslands, very short grazing lawns are present (Thapa et al., 2021). In the mixed tall grasslands (>2 m) the Narenga porphyrocoma and Erianthus ravennae species are abundant and these grasslands are considered to lie on a continuum with the short grasslands (Peet et al., 1999b).

3 Data and methods

3.1 Outline of approach

Figure 2 provides a flow chart of the methodology. We used topographic maps (1927, 1984) and aerial imagery of 1964 to assess the historic development of vegetation in the park and used LANDSAT imagery for reconstruction of annual land cover maps between 1993 and 2019. During the post-monsoon of 2019, BNP was visited to collect ground truth data of vegetation types. This dataset was used for supervised classification of land cover using a random forest model. From this time series, the areal extent, transitions and pattern of land cover were calculated. These results were related to environmental variables (precipitation, discharge, fire occurrences) to provide insight into the dynamics of land cover and the response to variation in environmental variables.

3.2 Satellite imagery and aerial photographs

To establish a time series of yearly land cover, we selected imagery of the Surface Reflectance dataset (Tier 1) of the LANDSAT 5, 7 and 8 satellites between 1993 and 2019. This imagery has a resolution of 30 m, and the full pre-processing of the images is elaborated in Appendix B. The composites of 1995, 1997 and 2006 were not classified. Too few images of sufficient quality were available for 1995 and 1997, and a too large number of erroneous pixel values in the image composite was observed for 2006, caused by the LANDSAT 7 Scan-line error. For each year we calculated Normalized Difference Vegetation Index (NDVI) (Rouse et al., 1974; Rock et al., 1986; Myneni et al., 1995; Beeri et al., 2007), Normalized Difference Drought Index (NDDI) (Gu et al., 2007; Biswas, 2010) and tasselled cap transformations (Crist and Cicone, 1984). We also included seasonal information by forming seasonal composites and a harmonic trend analysis of NDVI. In addition, aerial photographs from 1964 were obtained from the Forest Research and Training Centre, Ministry of Forests and Environment, Kathmandu. The individual photographs were combined, and a point cloud was generated, followed by an orthomosaic with a resolution of 2–3 m (Fig. S1). The resulting orthomosaic was resampled to a resolution of 30 m to enable comparison with the annual land cover time series.

Table 1 Datasets of the environmental drivers (a, b) and land cover dynamics (c), the derived metrics, and their time span.

Table 2 The metrics that are used for establishing effects of drivers on land cover dynamics, including the method and the part of the study area that is considered, and fragmentation metrics used for quantification of the grassland pattern (McGarigal and Marks, 1995; McGarigal et al., 2002; Plexida et al., 2014). The lower half of the table is partly adapted from Sertel et al. (2018).

3.3 Field data

Vegetation types were classified according to two schemes based on vegetation assemblages in BNP (Appendix C) (Dinerstein, 1979b; Jnawali and Wegge, 1993; Peet et al., 1999a) and remote sensing studies in similar environments (Dinerstein, 1979a; Sharma, 1999; Sarma et al., 2008; Biswas, 2010; Arieira et al., 2011; Thapa, 2011; Biswas et al., 2014). The first land cover series (level 1) contains four classes (water, bare substrate, grassland and forest) and spans from 1993 to 2019. The second series (level 2) uses solely LANDSAT 8 imagery (2013–2019). The eight classes consist of water, bare substrate, alluvial tall grasslands, short grasslands, mixed tall grassland, shrubland, riverine forest and Sal forest (Appendix C). The height and cover of vegetation, together with the dominant species present, were used to assign the vegetation types (Appendix D). These criteria enable discrimination between riverine grassland and non-riverine grasslands in the TAL (Biswas et al., 2014). For short grasslands the height criterion of <2 m is used. Shrublands (<5 m) entail no distinct vegetation assemblage recorded in BNP but were used in the TAL by Biswas (2010) and Sarma et al. (2008) and provided useful successional information on the impact of disturbances in Canada (Hermosilla et al., 2018).

3.4 Classification and validation

The random forest model (Breiman, 2001) was used for classification, considering its earlier successful use in mapping of vegetation in and near floodplains (Belgiu and Drãgu, 2016; Van Iersel et al., 2018; Harezlak et al., 2020). Two time series of land cover maps (level 1 and level 2) were created, and the aerial photographs from 1964 were classified manually (level 1) (Fig. S1). For 2019, the classification model was trained on the image composites of 2019 with 70 % of the field samples, and the other 30 % were used for validation. An error matrix was constructed for accuracy assessment and the user's and producer's accuracy was calculated (FAO, 2016). For classifying earlier years, all field samples were used to train a classification model on the two seasonal image composites of 2019 and image-derived attributes. The trained model for 2019 is transferable through time (Gómez et al., 2016). Pixels of each LANDSAT were then labelled with the most voted class for these pixels (Pal, 2005; Lawrence et al., 2006). To validate earlier years in the level 1 time series, additional validation samples (106) were collected with publicly available imagery via Google Earth for the available years (2000, 2010, 2011 and 2018). The imagery was not detailed enough for delineation of the level 2 classes. A grid was created digitally with a distance of 1200 m between sample locations outside the floodplain and 600 m inside the floodplain.

3.5 Environmental drivers

3.5.1 Meteorologic and hydrologic data

Hydrologic and meteorologic datasets of the Karnali River were obtained from the Department of Hydrology and Meteorology office in Kathmandu. The datasets include the yearly maximum and minimum discharges, daily discharges, and monthly precipitation. These data were measured at the Chisapani weather station north of the bifurcation and used to calculate hydrologic metrics (Table 1). Discharge has not been measured for the separate Karnali branches. Based on the dataset of 54 years of yearly peak discharges (Q_peak), extreme discharges (Q_ext) with ≥ 5-year recurrence times (≥ 12 500 m³ s⁻¹) were identified using the Gumbel distribution (Gumbel, 1958). It has been postulated that discharges with these recurrences cause major changes in the river courses in the Gangetic Plain (Thorne et al., 1993). The extreme precipitation events (P_ext) were identified by considering the peak rainfall months for each year (P_peak). The switch of the dominant discharge branch of the Karnali River system at the megafan was identified (Sinclair et al., 2017; Dingle et al., 2020b). The position of the water-filled channels during the dry period was extracted from the annual land cover maps (Table 1). We performed change detection of these water-filled channels by considering the land cover maps of time steps t and t+1, where t is 1 year. The results include newly covered area by water, disappeared water area and unchanged channels.

3.5.2 Forest fires

A forest fire dataset from 2000 to 2019 was obtained from NASA's Fire Information for Resource Management System (FIRMS). The dataset contains detections of fire derived from sensors on satellites (MODIS and VIIRS). Pixels that detected fire are flagged and attributes such as the brightness, confidence of occurrence and time of day are added. The number of pixels that were flagged with presence of a fire were selected for each calendar year, and we interpreted this number as the area that burned. This enables identification of years with high forest fire activity. It should be noted that this method results in an overestimation of the area that actually burned. The coarse resolutions (1 km for MODIS, 375 m for VIIRS) limit information on the exact extent of the fire. When a fire occurs at the border of two pixels, the two separate pixels are flagged, doubling its detected area. Also, when two fires occur at the same time on a spatial scale that is smaller than the resolution of the imagery, they are flagged as a single detection. The metric of forest fires was expressed as the percentage of area that is affected by fires for both the area within and outside the floodplain. For these two zones, the surface area is corrected for the areal extent of water and bare substrate.

3.6 Analysis of land cover dynamics

From the level 1 time series, two metrics were derived representing land cover dynamics (Table 2). Firstly, interannual transitions between land cover classes were calculated on a pixel basis with yearly time steps (Mas and Vega, 2012). In total, 16 combinations were calculated for the level 1 classification, with 64 combinations calculated for the level 2 classification. From these combinations we identified “removal of vegetation”, which is the transition from forest and grassland towards bare substrate and water, and the “transition of grassland to bare substrate” (Table 2). Secondly, the land cover dynamics were characterized by their pattern with landscape fragmentation metrics. These metrics enable quantification of spatial patterns through time (McGarigal and Marks, 1995; McGarigal et al., 2002; Plexida et al., 2014) and have been used for similar ecosystems in the Terai (Biswas, 2010; Thapa, 2011). We then identified years where extremes in environmental drivers (Q_ext and P_ext) coincide with large interannual transitions of grassland to bare substrate (Table 2). This was substantiated with statistical tests (Spearman tests in scatterplots) between the time series of annual extremes of environmental drivers (Q_peak, P_peak) and the interannual transitions of land cover (Table 2).

4 Results

The spatio-temporal dynamics of the 27-year land cover series are first presented alongside the historic analysis of land cover in 1927, 1964 and 1984. Hereafter, the datasets on environmental drivers are presented and coupled to the yearly spatial extent of land cover classes and transitions from grassland to bare substrate.

4.1 Spatio-temporal land cover dynamics in BNP

4.1.1 Land cover maps

For the reference year 2019, the overall accuracies are 84.7 % for the level 1 classification and 75.7 % for the level 2 classification (Appendix E). The classes of mixed tall grassland and shrubland underperform for the level 2 classification. For the level 1 classification, the calculated accuracies were 75.5 % (2000), 86.8 % (2010), 84.9 % (2011) and 78.3 % (2018) (Appendix E). Classes that underperform are bare substrate for the user's accuracy and water and grassland for the producer's accuracy.

A significant change in land cover (deforestation) around BNP can be seen when comparing the topographic maps of 1927 and 1984, particularly at the area between the Kauriala and Geruwa Rivers (Supplement S2). A number of grassland patches observed between 1993 and 2019 can be traced back to 1964 (Fig. 3f). However, grasslands have shrunk in size and surrounding smaller patches completely disappeared (Fig. 3f). The Baghaura, Lamkauli and Khauraha phantas are part of the remaining grasslands and consist of short and mixed tall grasslands (Fig. 3). The Lamkauli grassland can also be observed on the topographic map of 1927 as the village, settlement or grassland of Lathwa. No signs of a village or agricultural fields are observed on the 1964 imagery, in line with the record that the Baghaura and Lamkauli phantas were cultivated after 1965 (Pokheral, 1993; Brown, 1997).

When traversing the study area from west to east, the classes encountered most abundantly are, in turn, bare substrate, alluvial tall grasslands, short and mixed tall grasslands, shrubland, riverine forest, and lastly Sal forest (Fig. 3d and e). Shrubland is located close to the perennial streams and the national highway that crosses the park from Chisapani to the eastern border of the study area. Sal forest forms a sharp boundary with other vegetation types, coinciding with a difference in elevation that decreases southwards (Fig. A2). Alluvial tall grasslands are dominant in and close to the active river channels and ephemeral streams, whereas shorter grasslands and mixed tall grasslands are located more distal from the active channel belt (Fig. 3d and e). The active stream channels of the Geruwa River migrated westwards between 1964 and 2019 (Fig. 3a and b). Migration of individual channels is also observed, for example the change in the dominant discharge branch in the northwest (Fig. 3a, b and c) and the abandonment a river channel southwest of the Baghaura phanta at 28.49^∘ N, 81.24^∘ E. This channel transitioned from an active river channel in 1964 (Fig. 3a) to grassland in 2000 (Fig. 3b) and to forest in 2019 (Fig. 3c). Between 1993 and 2019, grasslands and low-flow river channels have covered practically the entire active channel belt at least once, highlighting the dynamic nature of the braided river system (Fig. 3g). This indicates that the larger part of the riverine forest located in and close to the stream channels is not older than 27 years. Compared to alluvial tall grassland along the Geruwa River and ephemeral streams, the more distal short and mixed tall grasslands are present for a larger time span (Fig. 3h).

Figure 3 Overview of land cover maps. (a) The manual delineation of land cover classes from the aerial photograph composite of 1964, (b, c) the level 1 dataset, (d, e) the level 2 dataset, (f) the grasslands that are present in 1964 and during 2017–2019 and where they appeared and disappeared, (g, h) and the number of years in total that water and grassland were present during the time period 1993–2019 (calculated from level 1 dataset). The full series of land cover maps is presented in Fig. A3.

Figure 4 The temporal changes in land cover classes and pattern from the (a) level 1 dataset and (b, c) level 2 dataset. White columns represent years with no data. Water can be interpreted as average water cover during the dry season. Note that a maximum of 25 % was set in (c) for the surface area outside the floodplain for visualization purposes; the remaining 75 % of the area consists solely of Sal forest. (d, e, f, g) Landscape fragmentation metrics (Table 2) through time for the grassland class in the floodplain, derived from the level 1 dataset.

Figure 5 Environmental drivers. (a) Yearly maximum peak discharge (Q_peak) and identified extremes (Q_ext) are shown for the discharge dataset. (b) For each year the month with the most precipitation (P_max) is shown, along with the identified extreme years (P_ext). (c, d) Surface area affected by fires in and outside the floodplain as detected from the MODIS (c) and VIIRS (d) datasets, respectively. (e) Change in surface area (hectares) covered by channels for each year compared to previous year, where newly covered means area covered by new channels, disappeared means area abandoned by channels, and unchanged means area covered by channels in the current year and previous year.

Figure 6 Environmental variables related to land cover dynamics derived from the level 1 classification series. (a, b) Annual transitions from grassland to bare substrate in hectares (orange bars) and years that experienced extreme discharge events (>12 500 m³ s⁻¹) and peak precipitation events within the floodplain (a) and outside the floodplain (b). Panels (c), (d) and (e) show correlation tests of environmental and land cover variables. Panels (c) and (d) show the annual removal of vegetation (the transition of Sal forest and grassland towards water and bare substrate) as related to the magnitude of yearly peak discharge before the channel relocation (c) and after the channel relocation (d). Panel (e) shows the yearly peak rainfall events (millimetres in a month) and the transition of grassland to bare substrate. Note the difference in scale on the y axes.

4.1.2 Land cover dynamics 1993–2019

Within the floodplain, succession and retrogression are observed in the level 1 classification (Fig. 4a) from which three periods can be distinguished: 1993–1999, 2001–2008 and 2010–2019. These periods are characterized by a gradual decline in the areal extent of bare substrate and an increase in vegetation, particularly grassland. The three periods are separated by years wherein the surface area of bare substrate increased considerably as compared to the year before (2000, 2008 and 2009). In the first period, grassland is in an increasing trend until 1998. In this period, water covers on average 18 % of the floodplain, amounting to 932 ha. In the second period (2001–2008), the area covered by water (low-flow channels) was slightly larger (1020 ha on average), while forest gradually decreased and grasslands saw an increase. In the third period (2010–2019), the area of bare substrate and water decreased to the lowest value since 1993. Bare substrate decreased gradually, whereas the coverage of water halved more abruptly to 478 ha on average (−53 %) as compared to the second period (2001–2008). The level 2 classification confirms the expansion of vegetation (Fig. 4b). This expansion is most evident in the floodplain where the surface area of alluvial tall grassland increased from 280 to 1185 ha. The short grasses and mixed tall grasslands do not show such a trend. A notable decline in forest is observed in 2015, both on the floodplain (Fig. 4b, to grassland) and outside of the floodplain (Fig. 4c). The grassland pattern shows a decrease in heterogeneity in the landscape: the number of patches decreased (patch density, PD, Fig. 4d), the length of the edges of grassland with other classes decreased (edge density, ED, Fig. 4e), the patches have become more compact (landscape shape index, LSI, Fig. 4f) and the patches have become more connected to each other (aggregation index, AI, Fig. 4g).

4.2 Environmental drivers

The discharge regime of the Karnali near Chisapani is characterized from 1993–2015 by a large seasonal variation, with a mean discharge of 1389 m³ s⁻¹, a minimum of 173 m³ s⁻¹ and a maximum of 21 700 m³ s⁻¹ in August 2014. The median is 678 m³ s⁻¹. Extreme discharges of the Karnali River can be up to 3 times the yearly average peak discharge. The monthly precipitation dataset shows a similar monsoonal variation, with a mean monthly rainfall of 196 mm and a maximum of 1673 mm, a minimum of 0 mm and a median of 40 mm. For the studied period of 1993–2019, Fig. 5a shows the extreme peak discharges (Q_ext) with a recurrence time of more than 5 years as calculated with the Gumbel distribution, and it also displays the periods with extreme precipitation (P_ext in Fig. 5b). An overview of the daily discharge measurements is shown in Fig. A4. The peak discharges of the daily and yearly dataset do not have the exact same values but are in statistical correspondence with each other. We based our analysis on the yearly discharge dataset. A set of 2 months with extreme conditions were identified in 2007 and 2014, of which the second-largest event, with 1560 mm for the monthly maximum (in 2007), is 47 % larger than the third-largest event (Fig. 5b). On an annual basis, the highest cumulative precipitation was recorded for the same years: in 2007 and 2014 with 3293 and 3390 mm, respectively.

Forest fires were most abundant in the study area in the years 2012, 2014, 2016 and 2019 based on both the MODIS and VIIRS detection datasets (Fig. 5c and d). As of 2012 and onwards, the variance in the yearly affected area by fires in the floodplain increased: the maximum surface area affected increased and alternated with years wherein the surface area affected was very low. The channel dynamics were quantified by comparing the annual location of the stream channels in the dry period. The newly covered surface area of channels has decreased considerably after 2009 as compared to the years before 2009 (Fig. 5e). Most changes in surface area took place in the years 2007, 2008 and 2009, whereafter the average absolute area that changed decreased.

4.3 Relating environmental drivers to land cover dynamics

For the years 2000, 2008 and 2009, when extreme peak discharges occurred, a response in vegetation is observed in the floodplain in the form of increased transitions from grassland to bare substrate (Fig. 6a). These increased transitions are not observed for 2013 and 2014, while in these years the recorded peak discharges were higher than in 2000, 2008 and 2009. The major change of the river course occurred in 2009, when the dominant discharge branch relocated from the Geruwa River (in our study area) to the Kauriala River. The extreme discharges in 2009 most likely have markedly changed the channel geometry, with a decrease in retrogression events after its relocation. This is also supported by the results of areal development of land cover (Fig. 4): before the channel relocation, the area of bare substrate only increases in the years that experienced extreme discharges and surprisingly not after the channel relocation in 2013 and 2014. Outside the floodplain (Fig. 6b) and close to ephemeral streams the large transitions of grassland to bare substrate in 2007 and 2014 coincide with years that experienced extreme precipitation events. These transitions occurred mostly near the ephemeral streams on a much smaller scale (20 times) as compared to the floodplain of the Karnali River.

For yearly peak discharges, a significant (p<0.05) positive correlation is present between the maximum discharge in a year and the areal change from grassland to bare substrate within the floodplain in the period between 1993 and 2009 (Fig. 6c). This is not the case for 2010–2019 (p>0.05), where the high peak discharges of 2013 and 2014 experience no increased transitions from grassland to bare substrate compared to years with average peak discharges (Fig. 6d). Outside of the floodplain, extreme rainfall events and grassland to bare substrate transitions show no correlation (p>0.05), but the extreme precipitation months of 2007 and 2014 coincide with the highest transitions from grassland to bare (labelled in Fig. 6e).

5 Discussion

We analysed the spatio-temporal dynamics of land cover in Bardia National Park and the effects of environmental drivers and variation therein on these dynamics. The results are discussed and compared to similar ecosystems at the foot of mountain ranges. As land cover acts as a biophysical indicator of the state of an ecosystem (Henderson-Sellers and Pitman, 1992; Meyer and Turner, 1992; Hansen and Loveland, 2012), the implications for BNP and possible cascading consequences for the ecosystem and nature management are also considered.

5.1 Evaluation of obtained accuracies and limitations

The observed temporal trends in land cover echo natural trajectories of succession, apart from the rapid forest recovery observed in the years 1999–2000 and the unrealistic fluctuations in riverine and Sal forest for the level 2 classification. The accuracies are similar to contemporary studies (Hermosilla et al., 2018; Harezlak et al., 2020; Van Iersel, 2020). Therefore, we argue that the series of land cover maps are of sufficient accuracy (76 %–87 % for level 1, 75 % for level 2) to draw conclusions on the general trends in land cover dynamics in the Geruwa River floodplain and adjacent part of Bardia National Park. The decreased surface area of grasslands is in accordance with observations in previous studies (Peet et al., 1999b; Jnawali and Wegge, 2000; Odden and Wegge, 2005; Neupane et al., 2020). Between 1976 and 1997, Odden (2007) recorded a decrease in grassland (−8 %) in the southwest of the study area. Between 1990 and 2013, Neupane et al. (2020) recorded a decrease in the surface area of grasslands from 31 % of the total surface area in BNP to 17 % (−45 %). In this study we recorded a decrease of −35 % between 1964 and 1996. The differences in numbers could be due to (1) the studies differing in terms of the time span studied, (2) the studies differing in terms of the area studied, (3) the studied years being in different phases of the biogeomorphologic cycle (Corenblit et al., 2007), causing the areal extent of floodplain grasslands to differ greatly prior to or after a flood (this cycle is also reflected in our results; Fig. 4), or (4) the methods differing in their approach to constructing land cover maps. Uncertainties in our results could originate from the classification models used, the input data (field samples and satellite image composites), validation data (collected in the field and with Google Earth), the change detection algorithm and the datasets of the environmental drivers.

5.2 Land cover dynamics from 1964 to 2019

The deforestation outside BNP is in line with the general observation that the TAL experienced high rates of deforestation in the previous century, but in contrast with the encroachment of forest within BNP in the past decades. The total areal extent of grassland was higher in 1964 (2614 ha) than in 2019 (2114 ha). Our interpretation is that short grasslands were proportionally more abundant in 1964 than in 2019 and vice versa for the alluvial tall grasslands, which expanded to 1184 ha in 2019 (Dinerstein, 1979b). Dinerstein (1979b) recorded that vegetation west and southwest of the Baghaura phanta was grassland with dispersed trees in 1976, which is in line with our findings for 1964. In 2019, this location is almost entirely covered by riverine forest (Fig. 3a and e). This is also the case for the Mansuri phanta, north of the Lamkauli phanta, which was a grassland up to 1976 but has converted to forest due to succession (Fig. 3f) (Bhatta, 2000).

The period between 1993 and 2019 shows gradual successional changes and years with abrupt retrogression events. A decrease in the heterogeneity of the grassland pattern is observed from 1993 to 2019 because grassland patches decreased in size and abundance and became more rounded and compacted. Riverine forest is observed in both the Karnali floodplain as well as along ephemeral streams. Bolton (1976) recorded the latter as an “undescribed riverine forest type”. Bare substrate is quickly colonized by early successional tall grassland. Concurrently, a portion of these grasslands transition towards the later successional grasslands and forests but at a slower rate. Lehmkuhl (1989) pointed out these different rates of succession for vegetation in Chitwan National Park. It is of interest to which extent the alluvial grasslands and riverine forest eventually expand, as in the last decade their combined area shows an increasing trend.

5.3 Environmental drivers

5.3.1 Hydrologic and meteorological drivers

The large discharge and precipitation fluctuations are associated with the monsoonal climate at the foot of mountain ranges (USAID, 2018) and alluvial megafans (Leier et al., 2005) such as the Karnali megafan and Kosi megafan. In 2014, the extreme precipitation in BNP and extreme discharge of the Karnali River coincided. Extensive loss of landmass and destruction of houses was recorded (Sinclair et al., 2017), and extreme flooding was also recorded in the nearby Babai River in BNP (Chhetri et al., 2020). During the peak discharge in 1983, second to the discharge in 2014 in magnitude, extreme flooding resulted in coverage of the Baghaura phanta with sediment (Bhatta, 2000). Between 1993 and 2019, we observed no conversion from grassland to bare substrate on this phanta. The active stream channels of the Geruwa River migrated westwards, which was also recorded by Rahhal et al. (2021).

For 2009 there are indications that a major redistribution of discharge occurred at the apex of the Karnali megafan. Dingle et al. (2020a) suggested, based on LANDSAT imagery, that the western branch was already dominant in 2001. This is not reflected in the results of the annual land cover time series of this study, which show a relocation of the main discharge in 2009. In addition, the rate of channel migration declined since 2009, and the water coverage in the Geruwa River decreased in the dry period (−53 %). In contrast, Rakhal et al. (2021) found an increase in long-term channel migration in the eastern branch of the Karnali River between 1977 and 2013. The cause for the discrepancy can possibly be found in the time span studied: the last decade (2010–2019) was not fully assessed by Rakhal et al. (2021) and the large-scale channel migrations in 2008 and 2009 could dominate the signal. During low discharges, the current distribution of discharge is considered to be about 80 % for the Kauriala River and 20 % for the Geruwa River (Sinclair et al., 2017; Dingle et al., 2020a). During peak flow of monsoonal discharges, the distribution is thought to be higher in the Kauriala River (55 %–65 %) (Dingle et al., 2017, 2020a). Timescales associated with avulsions at the apex of the Karnali megafan are estimated to be between ∼ 400 and 7000 years, (Dingle et al., 2020b), suggesting that the new discharge situation at the apex could hold for centuries. It is estimated that a discharge of 23 500 to 31 500 m³ s⁻¹ is required to entrain the coarse gravel near the apex of the Karnali megafan (Dingle et al., 2020b). This is higher than the discharge measured in 2009 (17 000 m³ s⁻¹) that caused the channel relocation. Altogether, indications are very strong that a major change in channel geometry, and consequently a redistribution of discharge, indeed occurred, resulting in more static hydromorphological processes in the Geruwa River. If erosion and sedimentation rates decrease, rejuvenation of the floodplain also decreases (Gunderson, 2002; Corenblit et al., 2007). The extreme peak discharges are not the sole hydrological driver for vegetation (Pickett et al., 1987; Corenblit et al., 2007). For example, prolonged inundation is common in a number of conservation areas in the TAL that are more characterized as wetland complexes, such as in Chitwan National Park (Dinerstein, 1979a). Our interpretation is that prolonged inundation is less of an important driver in BNP than in Chitwan National Park.

5.3.2 Forest fires, anthropogenic disturbances and herbivory

The variance of area affected by forest fires between years increased: the surface area affected (predominantly outside the floodplain) was large in the years 2012, 2014, 2016 and 2019, with increasing maxima. Ghimire et al. (2014) found that the surface area that burned increased after 2008, an increase that we observed since 2012. Ghimire et al. (2014) suggested changing climatic conditions as an explanation. To gain insight into the effect of anthropogenic activities, we interpreted the anthropogenic disturbances in relative terms, based on the days of access for locals, the yearly amount of permits issued for resource collection and what is known from literature on management activities (Appendix A) (Bolton, 1976; Brown, 1998; Bhatta, 2000; Bhattarai et al., 2017; Thapa et al., 2021). Anthropogenic disturbances by locals have strongly diminished since the establishment as a protected area in 1976, thereafter increasing towards the passing of the millennium, again followed by a considerable decrease (Bolton, 1976; Dinerstein, 1979b; Brown, 1997; Bhatta, 2000; Thapa et al., 2021). lease change to Persistent presence of (grazed) short grassland is still observed on the western bank of the Geruwa River, as confirmed by field observations on that location (Fig. A5). Maintenance by park staff also seem to have increased towards the passing of the millennium (Bhatta, 2000) and was extensive in recent years.

Herbivory affects land cover in a more gradual way than the pulse disturbances considered in this study, and therefore we do not attribute this driver to the abrupt vegetational changes. Over the long term, interpretation of grazing pressure on the vegetation pattern is hindered by uncertainty around the development of numbers of livestock and wild herbivores before and after the establishment of the conservation area. In 1975, biomass of livestock was 15–17 times higher (47 959 kg km⁻²) than biomass of wild ungulates in 1977 (2842–3120 kg km⁻²) in the southwestern part of the study area (Dinerstein, 1980). Between 1976 and 1999, grazing by livestock decreased, and ungulates are thought to have seen a 4-fold increase in numbers (Wegge and Storaas, 2009). It is of interest whether wild herbivores could have compensated for the decline in grazing pressure by livestock, on which no consensus exists in the literature (Karki et al., 2016; van Lunenburg et al., 2017; Wegge et al., 2019).

5.4 The effects of drivers on land cover dynamics

For the period before 2009, the results show a correlation between the transitions of grassland to bare substrate and the magnitude of yearly peak discharge. This relation is absent after 2009, although in 2013 and 2014 the peak discharges were of a higher magnitude than in 2009. This indicates that the impact of hydromorphological disturbances on vegetation in the floodplains of the Geruwa River decreased considerably during extreme discharges of the Karnali River. The absence of retrogression events during extreme discharges of the Karnali River in 2013 and 2014 contrasts with our hypothesis that in Bardia National Park extreme discharges of the Karnali River always cause a removal of vegetation, followed by a period of regrowth. The reduced hydromorphological dynamics favour transition to later successional stages of vegetation in the floodplain (Corenblit et al., 2007). This is in accordance with our observed expansions of alluvial tall grasslands and to a lesser extent that of riverine forest after 2009. The two most extreme precipitation events coincide with the largest transitions from grassland to bare substrate in the area outside of the floodplain along the ephemeral streams. Only during extreme precipitation events the magnitude of erosion and sedimentation processes seem large enough to convert vegetation to bare land because no correlation was found between the yearly maxima of precipitation and the retrogression of grassland. On the spatial and temporal scale of this study, no signal was detected between the transitions of land cover classes and years with extensive forest fires. We interpret the forest fires to not be a principal factor that converts forest towards grassland, which was also suggested by Lehmkuhl (2000). However, fires do play an important role in maintaining grassland (Peet et al., 1999b; Bhatta, 2000; Lehmkuhl, 2000).

Altogether, we attribute the decrease in surface area and abundance of grassland patches first to the reduced anthropogenic disturbances by locals since the establishment as protected area in 1976 and second to the westward migration of channels in the Geruwa River that is followed by succession to later successional stages in the abandoned channels of the Geruwa River. The contrasting recent expansion of alluvial tall grasslands is attributed to a decline in hydromorphological dynamics in the Geruwa floodplain. Phantas such as Lamkauli, Baghaura and Khauraha experienced significantly less encroachment than surrounding smaller patches. In addition, a number of grasslands in the south of the study area have become more distinct in the land cover maps in recent years. This is attributed to removal of trees in the grassland patches because of maintenance practices (Fig. A3). These practices predominantly seem to focus on these phantas and less so on the smaller surrounding patches. This is attributed to anthropogenic maintenance practices that are predominantly executed on these phantas, in contrast with the surrounding smaller grassland patches. If the grasslands would have been left unmanaged, encroachment of the forest would have been more severe.

5.5 Global comparison

Disturbance-dependent grasslands, especially those located in high rainfall regions, are threatened by indigenous forest expansion, which has been recorded in Australia, Africa, India, North America and South America (Puyravaud et al., 2003; Banfai and Bowman, 2006; Silva et al., 2008; Wigley et al., 2009). Grasslands in the Nepalese TAL can be added to this list. The redistribution of discharge through the two branches is thought to have had a profound impact on land cover. A comparison is made to parable conservation areas near lowland alluvial megafans. In the Pantanal (Brazil), extreme hydrologic conditions are explicitly analysed as a driving force in the pattern of vegetation communities (Arieira et al., 2011). The Taquari megafan in the Pantanal experienced a shift in environmental drivers due to an avulsion, increasing susceptibility to deforestation and fire (Louzada et al., 2020). At the foot of the Himalayas in northeastern India, extensive vegetational changes towards dryer vegetation types occurred after shifts in river courses in Jaldapara Wildlife Sanctuary (Biswas et al., 2014) and in Manas National Park (Sarma et al., 2008). In Jaldapara Wildlife Sanctuary, floods in 1968 caused a relocation of the Torsa River. This is regarded as one of the primary drivers that caused a transition of grassland towards woodland (Biswas et al., 2014). For the Manas National Park, Sarma et al. (2008) stated that after the relocation of the river, the change of grassland to forest could have been partly prevented by more active management of the grasslands

These analogous scenarios demonstrate the possible consequences for BNP, namely a change in the (relative) dominance of environmental drivers and possibly a shift towards dryer vegetation types and an increased susceptibility to droughts and fires. Without increased disturbances from another driver, for example anthropogenic disturbances, one could expect that a proportion of the now extensive alluvial tall grassland will transition to a more homogenic (riverine) forest, especially considering the decreasing trend in the heterogeneity of grassland in the landscape. This evokes questions as to the degree that humans should intervene in a (semi-)natural system. The importance of “human-dependent” habitat is also interesting in the light of how grasslands are described in the literature. Classically, grasslands are seen as a transitional land cover from a successional point of view (Clements, 1916), while recent studies highlight the importance to subdivide grasslands from a conservation perspective (Buisson et al., 2019; Bond and Parr, 2010) into old-growth (or ancient) grasslands (Bond, 2008; Veldman et al., 2015) and cultural (or derived) grasslands. The old-growth grasslands are considered to harbour a higher biodiversity and endemism when compared to the cultural grasslands. In contrast to this subdivision for conservation goals, subtropical grasslands in the TAL that are remnant of anthropogenic occupation and cultivation (the phantas) serve an important ecosystem role for tiger conservation goals in Bardia National Park by sustaining prey populations.

5.6 From the perspective of megafauna

Migrating channels at the apex of the Karnali megafan and associated redistribution of discharge are connected to the land cover distribution in the western part of BNP and with it the habitat of threatened megafaunal species. With the recent ongoing increase in grassland near the Geruwa River, the spatial extent of favourable habitat for species that depend on alluvial tall grasslands, such as greater one-horned rhinoceros and hog deer (Jnawali and Wegge, 2000; Odden et al., 2005; Pradhan et al., 2008), increased, but the natural mechanism that maintains these early successional tall grasses has reduced in strength. In Chitwan, habitat heterogeneity was positively correlated with the abundance and habitat preference of prey species of tigers (Bhattarai and Kindlmann, 2012). Moe and Wegge (1994) suggested that a high density of ungulates in BNP is possible due to the fine mosaic pattern of habitats in the Karnali floodplain. However, heterogeneity of grasslands in BNP shows a decreasing trend. This could be unfavourable for chital deer and tigers, although their numbers suggest otherwise. Tiger population showed the largest increase in BNP among the conservation areas in Nepal (DNPWC and DFSC, 2022), particularly due to conservational efforts.

5.7 Future research

In the next few decades, changes in the environmental drivers are projected due to (among other reasons) changing climatic conditions. Precipitation, temperature and mean annual discharge are projected to increase in the Karnali River basin, except for during the post-monsoon, when smaller amounts of precipitation and discharge could be expected (Dahal et al., 2020; Rakhal et al., 2021). The combination of Earth observation with time series of environmental drivers adds to the understanding of spatio-temporal dynamics of land cover and proves to be useful for monitoring with minimal monetary costs. With future technological advances in mind, the algorithm could be improved upon (e.g. Hermosilla et al., 2018; Van Iersel, 2020) to enable detailed analysis of vegetation dynamics on the community or species level, including invasive species (Bhatta et al., 2020). In addition, further research into the drivers and their changing relative dominance will aid conservational efforts and management decisions in BNP and similar nature reserves where disturbance-dependent grasslands are under threat from forest intrusion.

6 Conclusions

We aimed to quantify the spatio-temporal land cover change and establish relationships with the environmental drivers. This was done in Bardia National Park (BNP), located along the highly dynamic Karnali River at the foot of the Himalayas in a climate that is conducive to forest growth. We focussed on disturbance-dependent grasslands that are vulnerable to encroachment. Grasslands in the floodplain of the Geruwa branch of the Karnali River and the adjacent part of BNP were more abundant and were located more eastwards in 1964 than between 1993 and 2019. Between 1993 and 2019, grasslands showed a steady decline in heterogeneity, pinpointing the decrease in size and abundance of the patches. The more recent ongoing increase in surface area of grasslands can be attributed to expansion of alluvial tall grasslands in the floodplain. The maximum surface area affected by forest fires increased, but at the same time the affected surface area became more variable between years. Since 2009, the Geruwa branch of the Karnali River was partially abandoned. Before 2009, the magnitude of discharge events shows a correlation with transitions from grassland to bare substrate. This correlation is not present anymore after 2009, and no increased transitions are observed for 2013 and 2014, while the recorded peak discharges were of a higher magnitude in these years. In addition, the channel dynamics of the eastern branch of the Karnali River have decreased since 2009. Hence, we confirm our hypothesis that from 1964 to 2019 a shift occurred from grasslands to riverine and Sal forest. This shift is predominantly attributed to a westward shift of the active stream channels of the Geruwa River and the establishment of the conservation area in 1976 with its associated reduced anthropogenic disturbances by locals. The more recent expansion of alluvial tall grasslands is attributed to the redistribution of discharge in 2009 of the Karnali River. In 2019, land cover in the floodplain was more favourable for faunal species that depend on alluvial tall grasslands such as hog deer (Hyelaphus porcinus) and greater one-horned rhinoceros (Rhinoceros unicornis), although the mechanism to maintain these grasslands reduced in strength. For faunal species that benefit from a heterogenic vegetation pattern, land cover is thought to have become less favourable between 1964 and 2019 as the heterogeneity and surface area of the grassland patches decreased. Because hydromorphological processes in the floodplain have become more static, other agents of disturbances are more paramount to prevent encroachment of grasslands.

Appendix A: Supporting figures

Figure A1 View from the tower (Bagh Machan). As entering these grasslands is not recommended, this view is used to assign ground truth locations for short grassland (patch on right-hand side), alluvial tall grassland (in front, Saccharum spontaneum dominated), and mixed tall grasslands and riverine forest in the background. The location is northwest from the Khauraha phanta. Ground truthing was done in combination with satellite imagery to obtain the exact coordinates. The photograph was taken by the authors.

Figure A2 Photograph of an incision of the Karnali River. In the outcrop a loamy layer is positioned on top of a more sandy layer containing rounded pebbles, cobbles and boulders. Sal forest can be seen on top of the cliff. On the right-hand side of the image tall grasses grow amongst the grasses classified in this study as short grassland. The cliff marks the spatial extent of fluvial disturbances of the Karnali River. The difference in elevation is largest in the north of the study area and gradually decreases southwards along the interface of riverine and Sal forest. The photograph was taken by the authors.

Figure A3 Land cover maps of the level 1 and level 2 series.

Figure A4 Hydrological year showing the daily discharge measurements and their calendar day. Red dots represent the identified extreme discharges (Q_ext) that are higher than 12 500 m³ s⁻¹.

Figure A5 Grazing of cattle on the western bank of the Geruwa River at the southwestern border of study area. This image demonstrates the grazing pressure herbivores pose on the vegetation pattern at the border of the park both at present and in earlier times within the park when cattle grazing was allowed. The photograph was taken by the authors.

Table A1 Overview of anthropogenic activities in BNP as recorded in the literature.

Appendix B: Pre-processing of imagery

To establish an annual time series of land cover, we selected imagery of the Surface Reflectance dataset (Tier 1) of the LANDSAT 5, 7 and 8 satellites between 1993 and 2019. This dataset of imagery is corrected for variation in the energy source (e.g. sun angle) and atmospheric effects (e.g. aerosol scattering and thin clouds) and is therefore best suited to use for diachronic analysis (Zanter, 2019; Young et al., 2017). We accounted for the differences in sensors by using coefficients (Roy et al., 2016) to correct LANDSAT OLI imagery to LANDSAT 7 ETM+ imagery. The incorporation of seasonal information improves the classification results, especially if height data are absent, as the spectral signatures of the separate classes are more distinct due to phenological differences of vegetation throughout the year (Van Iersel et al., 2016; Kelley et al., 2018). For each year, two separate collections of images were selected based on the day of the year (between 240–366 of the post-monsoon and 1–150 of the subsequent pre-monsoon). In addition, images were selected based on cloud cover (<50 %) and quality (>9) and a geometric RMSE (root-mean-square error) of <10 m, and we performed cloud masking using the Fmask algorithm (Zhu et al., 2015). To enhance the available spectral information, NDVI (Rouse et al., 1974; Rock et al., 1986; Myneni et al., 1995; Beeri et al., 2007), NDDI (Gu et al., 2007) and tasselled cap transformations (Crist and Cicone, 1984) were calculated, which have been shown to improve classification results of vegetation (Price et al., 2002; Biswas, 2010). From each seasonal selection the images were combined into a composite image by taking the median value of each seasonal collection. In addition, for each year a separate annual collection of images was created, which was used to calculate the seasonal variation in the NDVI metric with a harmonic trend for each pixel in terms of magnitude and timing of the annual change. This was added as an indicator for the supervised classification model. To minimize noise caused by clouds, artefacts and the LANDSAT 7 scan-line error, we repeated the seasonal selection of imagery five times for each year with slightly different selection dates. After classification, the modus of the resulting five land cover maps was taken to obtain the most prevailing classification for each pixel for each year. The composites of 1995, 1997 and 2006 were not classified. Too many clouds were present for 1995 and 1997, and a too large number of erroneous pixel values in the image composite was observed for 2006, caused by the LANDSAT 7 scan-line error. We used the Google Earth Engine (Gorelick et al., 2017), earlier applied for mapping changes in floodplain ecosystems (Harezlak et al., 2020; Zurqani et al., 2018; Donchyts et al., 2016; Van Iersel, 2020), for selecting, processing and classifying LANDSAT imagery, whereas RStudio was used for post-classification analysis.

Appendix C: Grouping of vegetation types

Table C1 Grouping of vegetation types. Vegetation assemblages present in Bardia National Park and their associated vegetation classes as used in the level 2 classification model in this study. Grass assemblages are from Peet et al. (1999a). Forest assemblages are from Dinerstein (1979a) and a modification of data from Jnawali and Wegge (1993) and Pokheral (1993).

Appendix D: Criteria for classification

Table D1 Dominant species per vegetation type used for assigning classes. Cover and height are used for further discrimination of classes.

Appendix E: Error matrices

Table E1 Confusion matrices of the level 1 and level 2 classifications for 2019, as calculated with 30 % of the ground truth data and confusion matrices for 2000, 2010, 2011 and 2018 based on a separate validation set (106 samples).