1. Introduction
The global human population is projected to grow to 9.7 billion by 2050 (Kee et al., 2021). Increased crop production and changes in consumption patterns are prerequisite for meeting future food demands. Agricultural crop production needs to increase by around 70% by 2050 to meet the growing demand of population. Most of the people of Bangladesh are dependent on agricultural activities; agriculture is considered one of the largest producing sectors of the economy comprising 13.31% of the GDP (World Bank, 2022). The performance of agricultural sector exerts an overwhelming impact on major macroeconomic objectives like employment generation, poverty alleviation, food security, and other economic and social forces (Daboh & Jackson, 2023). The total area of Bangladesh is about 14.8 million hectares of which 7.83 million is the Net Cultivable Area (NCA) (World Bank, 2022). The percentage of cultivable lands has been decreasing in Bangladesh due to its increasing use for settlements, roads, industries, urban and other infrastructure development. About 25% of the total area is not available for agricultural practices. The national demand for food grain production is increasing every year against the alarmingly decreasing trend of NCA. This is one of the major challenges for Bangladesh, including the issues of national food security. This challenge can be overcome through judicious and scientific utilization of the country’s agricultural lands (Kee et al., 2021; Daboh & Jackson, 2023).
Several initiatives, including Good Agricultural Practices (GAPs), the 2030 Agenda for Sustainable Development, have been developed throughout the world to increase crop production (Beyer & Wacker, 2022). The proposed strategies were based on transitioning towards ecological cropping practices to substantially reduce external inputs (fertilizers and pesticides) and boost crop production without expanding cropland and further degrading natural resources (Lakatos, 2021). The Government of Bangladesh spent a significant amount to develop the overall agricultural sector keeping in view of the goals set out in the Seventh Five Year Plan and National Agriculture Policy (World Bank, 2022). Although Bangladesh Government took initiatives to boost agricultural crop production, it would be quite impossible to feed the increasing population without proper scientific management in next 20 years. For maximizing agricultural production and for the potential use of soil resources, a quick acquisition of related comprehensive information is necessary for scaling up processes. Use of Remote Sensing (RS) and Geographic Information System (GIS) can play active role in accomplishing these processes. However, RS and GIS are not automated decision-making systems, but these are tools to data acquisition which helps to generate relevant maps in support of the decision-making process (Kee et al., 2021; Lei, 2021). The first steps toward properly addressing soil and agriculture related issues are the spatial distribution of soil properties including metals, and RS together with GIS help to map the spatial distribution. Various geostatistical analysis methods allow analyzing spatial data and then predict the unsampled data’s location with high accuracy (Shankar & Monika, 2017; Han et al., 2021). The current study aims to assess soil contamination by metals using advanced GIS technology in the study areas located in six upazilas of Bangladesh.
2. Materials and Methods
2.1. Study Area
Six upazilas (e.g. Parbatipur, Ghatail, Tangail Sadar, Laksam, Mollahat and Bagher Para) from five districts were selected for the study based on soil contamination and agro-environmental criteria (such as major cropping pattern, land type and land use).
2.2. Soil Analysis
Five random soil locations were selected from each upazila to test the soil contamination degree. Thus, a total of thirty locations were selected. Collected soil samples were air-dried and ground to pass through a 2 mm-sieve and stored at approximately 4˚C in plastic bags until analysis (Sparks, 2003). The hydrometer method was used to determine particle size distribution (Page & Keeney, 1982). Soil pH was measured in soil suspension at 1:2.5 (soil: water) using the pH-meter. Soil electrical conductivity (EC) was measured in soil paste extract using an EC meter. The dichromate oxidation method was used to determine soil organic carbon (OC) according to the Walkley and Black procedure. Total nitrogen was measured by Micro-Kjeldahl’s method, and total phosphorus using spectrophotometer at 490 nm wavelength. Available phosphorus was determined by ascorbic acid blue color method (Olson et al., 2016). Available sulphur was determined turbidimetrically after extracting with calcium di-hydrogen phosphate (Page & Keeney, 1982). Soil samples were digested with nitric acid and hydrofluoric acid to determine the concentrations of metals (e.g. Cr, Co, Cu, Cd, Pb and Zn) (Sparks, 2003).
2.3. Development of Methodology
An internationally applicable approach was sought to optimize agricultural land utilization through the application of GIS for choosing appropriate cropping patterns by the farmers. Firstly, land/ crop suitability database was updated and validated. Social and economic information were incorporated with the existing crop data for preparing land suitability database of the major cropping patterns. A suitability assessment framework was developed including a GIS based multi-layers database followed by preparing the updated land suitability map based on the selected upazilas and major cropping patterns (Khan et al., 2018). Soil surveys were conducted to increase coverage in areas where the data were sparse or non-existent (Bartkowski et al., 2020; Heppenstall et al., 2021). Finally lands/ crop suitability assessment was done through multi-criteria analysis and matrix development (Naprstek & Smith, 2019).
2.4. Data Collection and Assemblage
Database development was done after collecting the necessary data, following the systematic approaches including 1) Database design, 2) Data capturing/digitization, 3) Data quality checking and populating, 4) Uploading database into GIS based land use suitability software (Naprstek & Smith, 2019). Database design was carried out considering the nature, format, source and quantity of data. Database consistency was checked following predefined methods described by researchers (Arostegui et al., 2006; Law & Collins, 2022). As part of the database design process, the Entity-Relationship (E-R) diagram was developed. For this, the entities were identified through the data analysis based on the existing database and attributes were set to each entity. Finally, relationship was developed among the entities (Law & Collins, 2022). Data entry and processing were done after designing the database. GIS data were captured and processed in ArcGIS; attribute data were processed in MS Access. The processed database was then incorporated into Geodatabase system of ArcGIS. Accuracy and consistency were ensured for data entry (Lei et al., 2016) by checking the location of farm, land type and other relevant information. Quality control of the data was done through graphical and visual interpretation. All attribute data and spatial data were processed using ArcGIS software. Major databases were developed such as Digital Elevation Model database, meteorology database, soil database, agriculture database, socioeconomic database, land and crop suitability database (Law & Collins, 2022; Lei et al., 2016).
2.5. Updating Land Type Map
Updating was done for Digital Elevation Model (DEM) and Land Type Map (LTM). DEM was updated using the interpolation method known as Inverse Distance Weighted (IDW). IDW was used to create the DEM with a spatial resolution of 100 m × 100 m. This estimated surface values for each cell using the value and distance of nearby points. The updated LTM was prepared using the updated DEM, Land and Soil Resources Utilization Guide and GIS. The mapping unit boundary did not have any physical demarcation of land type.
2.6. Remote Sensing and Image Processing
The supervised imaging classification of the study areas was carried out using the satellite image. The image preprocessing was done based on the radiometric and atmospheric calibrations using ENVI software (version 5.4). The satellite image was acquired to obtain a land type map, with a spatial resolution of 30 m. (Bonneu & Thomas-Agnan, 2009). The spatial distribution maps of land type, cropping pattern and DEM were carried out using the interpolation Kriging method (Bruno et al., 2010), which was involved in the geostatistical analyses in ArcGIS software 10.4. Descriptive statistical analyses were done to identify the normal distribution of the data using the Shapiro-Wilk test, and the upper and lower values and outlier values using SPSS (version 21). The unmeasured values of soil heavy metals were predicted using the semi-variogram models according to their accuracy.
2.7. Development of Model
Model was developed using satellite image processing, data analysis and ArcGIS (version 10.7) mapping for assessing soil contamination and crop suitability. This model was developed considering the multi-criteria analysis of the different aspects of edaphic, agro-climatic and economic factors (Chamling & Bera, 2020).
3. Results and Discussion
Satellite maps were produced for the study areas (Figure 1).
Figure 1. Satellite images of (a) Parbatipur, (b) Ghatail, (c) Tangail Sadar, (d) Laksam, (e) Mollahat and (f) Bagher Para upazila (Source: IRS PAN, 03 and Landsat 5 TM, 2019).
3.1. Soil Properties of Study Areas
Quantitative data were extracted using ArcGIS and satellite images for the study areas. Six upazilas had different soil textures. The soil textures of Parbatipur and Ghatail were clay whereas Tangail Sadar, Laksam, Mollahat and Bagher Para had silty loam textures. Soil pH ranged between 6.24 and 8.78, with an average of 712 ± 0.76 (Table 1). Regarding electrical conductivity (EC) values, descriptive data indicated that the study areas had a wide range of EC varying from 1.53 to 10.96 dS/m (6.47 ± 0.73 dS/m). The EC values < 2 dS/ m relatively occupied 16.6% of the study area and the rest 83.4% recorded EC values > 2 dS/m. According to the agronomic classification of soil salinity (Hammam & Mohamed, 2020), only one sample had EC < 2 dS/m (non-saline) and 95% of the samples classified as slightly saline soils (2 - 4 dS/m). The obtained data revealed that soil organic carbon varied between 4.34 and 18.92 (10.69 ± 1.85) mg/kg (Table 1). Maximum level of total nitrogen (34.34 ppm) was found in the soils of Laksam upazila. However, minimum content of total N (8.11 ppm) and total phosphorus (19.16 ppm) was present in the soils of Tangail Sadar. The content of total P was found to be the maximum in the soils of Parbatipur (38.64 ppm). Regarding the available forms of nutrients, maximum and minimum levels of P were found in Bagher Para (23.97 ppm) and Tangail Sadar (9.34 ppm) respectively. Available S content was minimum in Mollahat (0.73 ppm) and maximum in Tangail Sadar (5.18 ppm) among the study areas (Table 1).
Table 1. Descriptive statistics of the studied variables.
Upazila |
pH |
EC (dS/m) |
OC (g/kg) |
Total N (ppm) |
Total P (ppm) |
Available P (ppm) |
Available S (ppm) |
Min |
Max |
Min |
Max |
Min |
Max |
Min |
Max |
Min |
Max |
Min |
Max |
Min |
Max |
Parbatipur |
7.71 |
8.78 |
2.57 |
10.54 |
7.32 |
17.09 |
12.01 |
32.34 |
24.86 |
38.64 |
10.45 |
18.67 |
0.93 |
4.87 |
Ghatail |
6.24 |
7.13 |
1.53 |
9.75 |
4.87 |
15.85 |
10.89 |
33.23 |
20.86 |
32.64 |
11.35 |
19.07 |
1.04 |
3.87 |
Tangail Sadar |
6.89 |
7.26 |
2.65 |
10.96 |
4.34 |
14.75 |
8.11 |
29.03 |
19.16 |
29.64 |
9.34 |
17.07 |
0.94 |
5.18 |
Laksam |
5.98 |
6.24 |
2.65 |
10.13 |
6.38 |
18.92 |
10.35 |
34.34 |
24.34 |
31.42 |
10.15 |
20.17 |
0.79 |
2.62 |
Mollahat |
6.39 |
7.37 |
2.22 |
9.39 |
7.02 |
16.72 |
14.72 |
30.13 |
25.06 |
30.04 |
10.05 |
17.77 |
0.73 |
1.07 |
Bagher Para |
7.31 |
8.15 |
2.77 |
10.11 |
5.47 |
18.46 |
10.05 |
31.93 |
19.54 |
29.78 |
12.45 |
23.97 |
0.98 |
3.43 |
3.2. Metal Contents
Descriptive statistics for total concentrations of the investigated Cr, Co, Cu, Cd, Pb and Zn are detailed in Table 2. The Cr total concentration ranged between 36.31 and 164.06 with a mean of 101.70 ± 58.36 mg/kg. Cobalt concentration varied from 18.34 to 46.98 with an average of 32.65 ± 11.34 mg/kg. The range of total Cu concentration was between 24.31 and 91.24, with an average of 58.64 ± 10.07 mg/kg. The total content of Cd has an average of 1.25 ± 0.09 mg/kg. Total concentrations of Pb and Zn as microelements were 12.05 to 81.13 (46.10 ± 25.76) mg/kg, and 54.39 to 187.02 (104.65 ± 45.74) mg/kg, respectively. Statistical analysis observed that metal concentrations varied significantly across all study areas. Chromium concentrations in Tangail Sadar soils significantly (p ≤ 0.05) exceeded those in Parbatipur, Ghatail, Laksam, Mollahat and Bagher Para. The upazilas with significantly (p ≤ 0.05) highest Co and Cu concentrations were Parbatipur and Mollahat respectively. The Cd levels in Bagher Para soils were significantly (p ≤ 0.05) higher than in the Parbatipur, Ghatail, Laksam and Mollahat upazila soils. Laksam soils had significantly (p ≤ 0.05) highest Pb and Zn concentrations (Table 2).
Table 2. Metal concentrations in the soils of study areas.
Upazila |
Cr (mg/kg) |
Co (mg/kg) |
Cu (mg/kg) |
Cd (mg/kg) |
Pb (mg/kg) |
Zn (mg/kg) |
Min |
Max |
Min |
Max |
Min |
Max |
Min |
Max |
Min |
Max |
Min |
Max |
Parbatipur |
50.31 |
162.74 |
21.54 |
46.98 |
34.39 |
90.75 |
0.32 |
2.09 |
15.24 |
78.34 |
66.86 |
136.61 |
Ghatail |
42.11 |
156.74 |
20.57 |
41.08 |
30.93 |
89.54 |
0.26 |
1.59 |
20.72 |
72.13 |
60.96 |
124.34 |
Tangail Sadar |
51.31 |
164.06 |
18.34 |
37.93 |
30.05 |
87.97 |
0.34 |
2.35 |
19.05 |
65.93 |
58.13 |
149.38 |
Laksam |
36.31 |
140.24 |
27.73 |
45.68 |
24.31 |
78.15 |
0.24 |
2.56 |
19.72 |
81.13 |
54.39 |
187.02 |
Mollahat |
45.21 |
158.76 |
21.47 |
43.07 |
25.98 |
91.24 |
0.31 |
1.82 |
12.05 |
73.93 |
73.11 |
135.48 |
Bagher Para |
52.31 |
160.35 |
23.58 |
43.91 |
35.34 |
85.07 |
0.35 |
2.74 |
14.05 |
80.93 |
67.03 |
142.44 |
All means of the studied metals were higher than the chemical composition of the upper continental crust ranging between 0.12 and 94.65 which was observed by Bradl (2005). The wide range of metal concentrations in the study areas could be due to the accumulation of metals and metalloids through emissions from the rapidly expanding industrial areas, mine tailings, disposal of wastes, land application of chemical fertilizers, sewage sludge, pesticides, wastewater irrigation, coal combustion residues and spillage of petrochemicals (Khan et al., 2014). The obtained results agreed with other investigations in the previous literatures (Shelbaya et al., 2021; Hammam & Mohamed, 2020). Values of skewness (Pearson) for Cr ranged between 1.15 and 1.87 among the study areas. Likewise, values of skewness (Pearson) for Co, Cu, Cd, Pb and Zn were 0.42 - 0.67, 1.32 - 1.89, 0.81 - 1.03, 1.65 - 1.96, and 2.42 - 2.94, respectively (Table 3).
Table 3. Descriptive statistics of metals in the soils of study areas.
Upazila |
Cr (mg/kg) |
Co (mg/kg) |
Cu (mg/kg) |
Cd (mg/kg) |
Pb (mg/kg) |
Zn (mg/kg) |
Cr (%) |
Co (%) |
Cu (%) |
Cd (%) |
Pb (%) |
Zn (%) |
Skewness (Pearson) |
coefficient of variation (CV) |
Parbatipur |
1.03 |
0.47 |
1.43 |
0.89 |
1.65 |
2.78 |
82.02 |
42.76 |
98.27 |
16.95 |
60.23 |
51.87 |
Ghatail |
1.15 |
0.45 |
1.89 |
0.93 |
1.68 |
2.42 |
86.04 |
44.34 |
103.71 |
24.64 |
57.12 |
48.67 |
Tangail Sadar |
1.76 |
0.42 |
1.36 |
1.01 |
1.74 |
2.93 |
71.54 |
34.82 |
102.35 |
28.07 |
64.68 |
52.86 |
Laksam |
1.23 |
0.67 |
1.32 |
0.81 |
1.96 |
2.94 |
82.86 |
47.55 |
92.56 |
19.19 |
58.46 |
50.61 |
Mollahat |
1.87 |
0.59 |
1.72 |
0.92 |
1.82 |
2.91 |
79.12 |
35.08 |
100.86 |
18.75 |
62.18 |
45.09 |
Bagher Para |
1.42 |
0.54 |
1.45 |
1.03 |
1.94 |
2.51 |
74.17 |
41.59 |
102.32 |
26.05 |
63.02 |
54.81 |
To estimate the degree of influence from anthropogenic activities on soil metal concentrations, the coefficient of variation (CV) of a given metal in the soils has been used in many studies (Hammam & Mohamed, 2020). A high CV (>50%) suggested that there were larger variations in metal concentrations among the different sites and that spatial distribution of the metal concentrations in the study areas were non-homogeneous (Gao et al., 2010). The CVs of the metals in the soils of six upazilas were as follows: Cu (100.01% ± 4.09%) > Cr (79.29% ± 5.51%) > Pb (60.95% ± 2.87%) > Zn (50.65% ± 3.42%) > Co (41.02% ± 5.11%) > Cd (22.28% ± 4.55%) (Table 3). The CVs of metals indicate the distribution uniformity and degree of variation of metal elements in the study area (Shelbaya et al., 2021). The larger values of CV for Cu, Cr, Pb and Zn indicated greater influence of human activities (Gao et al., 2010). These CV values also suggested that the content of metals in soil was not only affected by the geological background but also by human activities (Hu et al., 2014).
3.3. Spatial Data Modeling
Present study indicated the parameters of semi-variogram modeling. The accuracy parameter of the best fit model among the semi-variograms was calculated using RMSE, MSE, RMSSE, and ASE. The results showed that the Gaussiang model was highly fit for soil pH and Pb; the Exponential model was fit for EC; the Stable model was fit for OC, Co, Cu, and Cd. Furthermore, the Spherical model was fit for both Cr and Zn. The MSE values were close to zero for all selected parameters, while the values of RMSE and RMSSE were close to one. However, values of ASE were small in all the resulted models.
3.4. Spatial Distribution of Cropping Pattern
The net cultivable area (NCA) and cropping pattern were determined using ArcGIS and ENVI software. Cropping pattern mapping was produced through crop monitoring, which can be based on crop type identification using multi-temporal data (Leblois et al., 2014). After collecting necessary data, updated cropping type maps of the selected upazilas (Figure 2) were validated. Ground truthing of the updated maps were found to be valid in all cases. In Tangail upazila, no cultivable land was found because of highly developed settlement. It was apparent from the study that the updated land which was cultivable only a year back was converted to a settlement area by a housing company. In Mollahat upazila, no valid highland was found. This discrepancy was due to land levelling for cultivation of shrimp in that area. Another selected study location in Parbatipur upazila, near the Maddhyapara Coal Mine, was adversely affected due to coal mining. The land was subsided significantly and consequently the medium lowland was turned into a perennial water body.
The NCA in Parbatipur upazila was found to be ~30,060 ha. Single, double and triple cropping areas were ~5%, 76% and 19% of the NCA respectively. The main cropping pattern was found Fallow-T. Aman-Boro covering ~69.1% of the NCA. The NCA of Ghatail upazila was found ~30,150 ha where double cropped area was dominant (~68%). About 18% of the NCA was covered by single cropping whereas the remaining 14% was covered by triple cropped lands. The maximum area was covered by Fallow-T. Aman-Boro cropping pattern that covered ~59% of the NCA. The second highest coverage was by perennial crops covering ~12% of the NCA. Among the perennial crops, banana, pineapple, ginger and turmeric covered ~1500, 1000, 150 and 950 ha respectively. The NCA of Tangail Sadar was found ~20,465 ha. Single, double and triple cropped patterns practiced were ~16%, 55% and 30% of the NCA respectively. The main cropping pattern was Fallow-Fallow-Rabi crops that occupied ~20.1% of the NCA. The other major cropping patterns were Fallow-Fallow-Boro, Fallow-T. Aman-Rabi crops, Fallow-T. Aman-Boro, and Jute-T. Aman-Boro covering ~12.4%, 11.6%, 11.2% and 11.6% of the NCA respectively. A small area was covered with single perennial crops like banana (~150 ha), papaya (~45 ha) and turmeric (~42 ha). In Laksam upazila, the NCA was ~11,010 ha, where triple cropping pattern was dominant. The dominant cropping pattern was observed Fallow-Aman-Boro which covered ~41.6% of the total NCA. The NCA of Mollahat upazila was ~12,130 ha. Single, double and triple cropped patterns were ~38%, 45% and 16% of the NCA respectively. The main cropping pattern was Fallow-Fallow-Boro which covered ~31.7% of the NCA. Among the Rabi crops, the percentage of coverage of pulses, vegetables, mustard, spices and wheat were ~10%, 7%, 6%, 1% and 1% of the NCA respectively. The NCA was ~19,020 ha in Bagher Para. Double and triple cropped areas were equally dominant in this upazila. The maximum area was covered by Fallow-T. Aman-Rabi crops/ Boro indicating ~50.1% of the NCA.
![]()
Figure 2. Updated cropping pattern map (Source: IRS PAN, 03 and Landsat 5 TM, 2019).
Mapping for crop type and cropping pattern has been studied in literatures (Woodward et al., 2012) using multi-temporal optical imagery. In the present study, we determined properties of soil (Table 1) and metal contents (Table 2) that ultimately impacted on crop type and cropping pattern. Our study is in agreement with previous studies (Gao et al., 2010; Hu et al., 2014; Woodward et al., 2012). Researchers determined biochemical, morphological, structural, physiological and phenological properties of soil influencing crop status, growth, reproduction and survival. Another group of researchers (Lillesand & Kiefer, 2000) considered crop calendar, which is vital in identifying different cropping patterns. The crop calendar information illustrates different cropping patterns, crop growing cycles at different times (Woodward et al., 2012). Researchers (Lillesand & Kiefer, 2000; El-Bastawesy et al., 2012) observed spatial distribution of different cropping patterns using remote sensing data. They used ERDAS Imagine software in remote sensing which functioned differently ENVI software used in our study. The functions of ERDAS such as the format conversion and image enhancement possess special tool boxes, and can parallelly operate, while the functions in ENVI are scattered in the form of a menu item, and are not able to undertake parallel processing (Dawood & Koshak, 2011).
4. Conclusion
The current study highlights the assessment of metal contamination in soil, which is considered one of the most important drawbacks to sustainable agricultural development. Cropping pattern is also a major factor contributing to sustainable development. The study suggested that metal concentrations varied significantly (p ≤ 0.05) across all study areas. The coefficient of variation (CV) for the selected metals had an order of Cu > Cr > Pb > Zn > Co > Cd, suggesting that the metals in soil were not only affected by the geological background but also by human activities. All means of the metals were higher than the chemical composition of the upper continental crust ranging between 0.12 and 94.65. The accuracy parameter of the best fit model among the semi-variograms was calculated using RMSE, MSE, RMSSE, and ASE. In Parbatipur, the main cropping pattern was found to be Fallow—T. Aman—Boro covering ~69.1% of the net cultivable area (NCA). Double cropped area was dominant (~68%) in Ghatail. In Mollahat, single, double and triple cropping patterns were ~38%, 45% and 16% of the NCA respectively. Double and triple cropped areas were equally dominant in Bagher Para.
Acknowledgements
The research was funded by the Ministry of Water Resources in Bangladesh. The authors would like to acknowledge Bangladesh Council of Scientific and Industrial Research (BCSIR) for providing necessary technical facilities. The authors are obliged to the Center for Environmental and Geographic Information Services (CEGIS) for providing necessary maps and satellite images.