2.1 Study site description
We studied the riparian zone along the Wang River in Thimphu City, Bhutan’s capital, and the largest city (Fig. 2a and b). It was the most rapidly developing city in Bhutan, with a spatial expansion rate of 36.3% between 2001 and 2017 (Yangzom et al., 2017). The city was home to 0.14 million of Bhutan's 0.7 million people in 2017, making it the most densely populated urban area in the country (National Statistics Bureau, 2018). The city is situated within a valley in the Himalayan Mountain range, constricted and narrowed by surrounding hills of blue pine forest (Pinus wallichiana A.B. Jacks.). Its topography is primarily gentle slopes along the riverbanks, gradually giving rise to steep slopes towards peri-urban zones. The Köppen-Geiger climate classification system classifies Thimphu as a subtropical highland climate, Cfb (Beck et al., 2018), characterised by distinct seasons with cold-dry winters and warm-wet summers (Vilà-Vilardell et al., 2020). The mean annual temperature is 14.7oC and the annual precipitation is 608.9 mm with peak seasonal rainfall (June to August) reaching 123.80 mm per month on average (Weather and Climate Services Division of Bhutan, 2018).
2.2 Sampling and sample collection
Sampling was conducted within 18 km of the RZ of the Wang River within the city’s urban boundary (Fig. 2c) in May 2023, before the onset of the monsoon. The RZ was defined as the part of the river channel between low and high water marks, hydrologically connected to the adjacent terrestrial landscape, and containing vegetation influenced by the flood regime (Naiman et al., 2010; Naiman et al., 1998). The narrow RZ in some sections of the city was completely replaced (Fig. 3a) or truncated (Fig. 3b) by IC, commonly comprising roads, flood protection walls, and concrete dikes. The lack of periodic flooding in the river valley has encouraged land holdings and habitation close to riverbanks, further narrowing the RZ.
The first site was located in Babesa at the southern end of the urban boundary. Subsequently, plots were established at intervals of ~600 m along the river until Kabesa, at the northern end of the urban boundary (Fig. 2c). Plots that fell within areas restricted to the public, such as the Royal Palace, were excluded. Each plot was field verified to ensure that data collection was feasible. Plots that were difficult or unsafe to access owing to steep terrain and fenced private property were relocated to the nearest feasible location. Sites with ongoing landscape development were also excluded, leaving a total of 25 sampling sites.
Sampling plots (n = 25) were grouped into two land-use classes, namely ‘urban land use’ (ULU) and ‘other land use’ (OLU), based on a land use land cover (LULC) map of the city (Fig. 4a). The LULC map was obtained by supervised classification of Sentinel-2 satellite images in the Google Earth Engine (GEE) using a random forest classifier (Truong et al., 2022). ULU primarily comprises ICs, such as buildings, roads, parking lots, recreational areas, and sidewalks. Vegetation, agriculture, and others were merged to form the OLU, comprising mainly parks, paddy fields, orchards, scrubs, and forests. To qualify plots as ULU or OLU, we used a 100 m radius circular buffer (Fig. 4c) in QGIS (3.28) to visually estimate land cover types from the LULC map of the city. The 100 m radius buffer was used because it covered a sufficient area to distinguish between different land use features, including the river channel (Fig. 4c). We then assigned plots with more than 20% IC coverage (using the IC coverage threshold adapted from Tippler et al. (2012)) to ULU sites, and the remainder to OLU sites. Historical images from Google Earth Pro (GEP) were used to verify site groupings, and similar visual estimations were observed, yielding ULU = 13 and OLU = 12 (Fig. 4).
At each site, a 5 × 5 m quadrat was established at a minimum distance of 1 m from the high-water mark, within which five samples of surface soil (0 – 20 cm depth, (Zhou et al., 2019)) were collected using an Edelman soil auger. Samples were taken from the quadrat’s four corners as well as the centre (Zhou et al., 2022). The samples from each site were homogenised (mixed) and stored in aluminium containers. The auger was cleaned with ethanol before obtaining samples from the new plots. Surface leaf litter was cleared before sample collection. A subsample of approximately 1 kg of soil collected from each site was used to analyse texture (as described by the National Committee on Soil and Terrain (2009)), soil pH (1:5 soil/water suspension test), total Kjeldahl nitrogen (TKN, semi-micro Kjeldahl), available phosphorus (P, Bray 1-P) and available potassium (K, CaCl2 extraction method) at the Soil and Plant Analytical Laboratory, National Soil Service Centre at Simtokha in Thimphu following methods from Rayment and Lyons (2011) and Houba et al. (2000). Another subsample (also approximately 1 kg) was used to analyse MP content (see s2.3).
The distance between the nearest IC and the sampling point was measured using a 100 m tape. The distance between each sampling point and the nearest upstream storm drain, wastewater or sewage discharge point was calculated using GEP. The percentage cover of surface plastic and textile waste in the 5 m × 5 m quadrat was visually estimated before sample collection.
2.3 MP extraction and quantification
The MP extraction method was adapted from Nematollahi et al. (2022). In the laboratory, each homogenised soil sample was placed in an aluminium container sealed with perforated aluminium foil to prevent cross-contamination and was oven-dried at 40 °C. Dried samples were sieved (stainless steel - 5 mm mesh), and any remaining larger rocks or plant debris were manually removed. 10 g of the sieved samples were mixed with 70 ml of 30% hydrogen peroxide (H2O2) in a 500 ml conical flask to digest the organic matter. Seven ml of H2O2 per gram of dried soil adequately covered the samples and allowed for an oxidation reaction. The samples were then placed in an oven at 40 °C to complete the digestion and dry any remaining H2O2.
The extraction of MP involved density separation with ZnCl2 solution (density = 1.6 gcm-3), which can separate most MP, including high-density plastics, from the soil matrix (Junhao et al., 2021). For each sample, 50 ml of ZnCl2 solution (filtered with a 2.0 µm Millipore glass fibre filter) was added and then shaken with an IKA digital orbital shaker (KS 501) at 250 rpm for 5 min (Li et al., 2020). After 24 h, the supernatant was transferred to a separate tube, centrifuged at 5000 rpm for 15 min, and vacuum filtered through a 2.0 µm Millipore glass fibre filter, while ZnCl2 was rinsed away with doubly distilled water (Nabi & Zhang, 2022). The filter was transferred to a Petri dish, covered with perforated aluminium foil, and oven-dried at 40 °C for 24 h (Yang et al., 2021). All soil residues were transferred to separate Petri dishes, covered with perforated aluminium foil, and oven-dried at 40 °C to check for the presence of MPs. Throughout the extraction process, the quality assurance recommendations of Dehghani et al. (2017), including blank sample preparation, Dehghani et al. (2017) were followed.
All filters and soil residues were observed under a stereo microscope (Zeiss SteREO Discovery.V20) equipped with a camera and computer connection. Each sample was checked twice for the presence of MP (Corradini et al., 2021). To verify MP and prevent the erroneous identification of soil mineral particles (such as mica) as MP (Fig. 5a and b), heat treatment was implemented (Y. Zhang et al., 2020). All MPs were detected on filter papers, and no MPs were found in soil residues or blanks. MPs were then counted and recorded by colour (Fig. 5c and d). The MPs obtained from each 10 g sample were extrapolated to a per kg scale to estimate the concentration of MPs in larger volumes.
2.4 Data Analysis
In this study, logistic regression analysis was employed to investigate the impact of land use, plot IC % and distance to IC on the probability of detecting MPs in RS. The MP count was transformed into a binary variable of presence/absence for this analysis. Additionally, a Mann-Whitney U test was conducted to compare differences in MP distribution, soil physiochemical properties, and anthropogenic variables between sites of ULU and OLU. Fisher’s exact test was utilized to explore the potential influence of soil texture on the presence of MPs in the RS. This test was deemed appropriate due to the small sample size, the presence of four categories of soil texture, and the binary response variable of MP particles.
Furthermore, Spearman’s correlation analysis was performed to assess the relationships between the abundance of MPs and various factors including soil pH, TKN, available P, available K, surface plastic waste, surface textile waste, and the distance to upstream drainage outlets. The Shapiro-Wilk test was used to evaluate the normality of the data. It was found that, except for the distance to the IC (p = 0.52), all other variables did not meet the assumption of normality (p < 0.05). Consequently, non-parametric statistical tests were employed for the analysis. All statistical analyses were conducted using IBM SPSS Statistics version 29. Maps depicting the spatial distribution of the study plots and land use classification were generated using QGIS version 3.28 and GEE, while BioRender was utilised to create graphical illustrations (https://app.biorender.com).