A New Multibranch Model for Metals in River Systems: Impacts and Control of Tannery Wastes in Bangladesh

: A new multibranch Integrated Catchment (INCA) model INCA-Metals has been developed to simulate the impact of tannery discharges on river systems. The model accounts for the key chemical reaction kinetic processes operating as well as sedimentation, resuspension, dilution, mixing and redistribution of pollutants in rivers downstream of tannery discharge points and for mine discharges or acid rock drainage sites. The model is dynamic and simulates the daily behaviour of hydrology and eight metals, including cadmium, mercury, copper, zinc, lead, arsenic, manganese and chromium, as well as cyanide and ammonia. The model is semi-distributed and can simulate catchments, tributaries and instream river behaviour. The model can also account for diffuse pollution from rural runoff as well as point sources from efﬂuent and trade discharges. The model has been applied to the new Savar tannery complex on the Dhaleshwari River system in Bangladesh to assess the impacts on pollution levels in the river system and to evaluate a set of treatment scenarios for pollution control, particularly in the dry season. It is shown that the new efﬂuent treatment plant at Savar needs to signiﬁcantly improve its operation and treatment capability in order to alleviate metal pollution in the downstream Dhaleshwari River System and also protect the Meghna River System that falls in the Bay of Bengal.


Introduction
Globally, polluted surface and ground water causes harm to ecosystems and human health. Sustainable Development Goal (SDG) 6 [1] for water and sanitation builds on earlier frameworks established by the Millennium Development Goals to address the complexity of interconnected water quality, quantity and affordability concerns. SDG indicator 6.3 seeks to achieve improvement in ambient water quality by significantly reducing flows of pollution into water bodies. In many developing countries where progress toward SDG 6 is off-target, combinations of untreated industrial, domestic, urban and agricultural waste contribute to dangerous levels of toxicity. The reduction of heavy metal pollution from leather tanneries, the garment industry and mining activities is a key issue and is of major concern for the protection of river water quality and ecology, as well as for sustaining livelihoods for people living close to river systems.
Metal pollution from tanneries is a major issue in countries such as Bangladesh and India [2,3], where metals such as chromium, copper and cadmium are released from the tannery processes, often with minimal effluent treatment, causing severe pollution of local streams and rivers. Other metal pollution arises from different industries such as mine effluents and acid rock drainage (ARD) which can have significant impacts on river systems [4]. Treatment methods for both tanneries and the mining industry are often problematic, as it is expensive, technically demanding and insufficiently monitored by regulatory authorities. Many active and passive waste treatment systems have been investigated and applied to control tannery effluents and mine discharges [5,6]. The problem of pollution following mine closure is particularly acute in many European countries where historical mining has been undertaken with poor pollution control and where the high costs of treatment make restoration a major issue [7].
This paper focuses on pollution in the Dhaleshwari River system in Dhaka, Bangladesh, where a new complex of tanneries at Savar was established in 2016 ( Figure 1). The Dhaleshwari river system is part of the Jamuna-Brahmaputra catchment and is the major source of freshwater to the Greater Dhaka Watershed. Flows in the river have decreased due to natural as well as human interventions, but augmentation of flow during the dry season in the Dhaleshwari and some other rivers is under consideration targeting pollution reduction and ecological restoration works planned under the Bangladesh River Masterplan and other initiatives. Metal pollution from tanneries is a major issue in countries such as Bangladesh and India [2,3], where metals such as chromium, copper and cadmium are released from the tannery processes, often with minimal effluent treatment, causing severe pollution of local streams and rivers. Other metal pollution arises from different industries such as mine effluents and acid rock drainage (ARD) which can have significant impacts on river systems [4]. Treatment methods for both tanneries and the mining industry are often problematic, as it is expensive, technically demanding and insufficiently monitored by regulatory authorities. Many active and passive waste treatment systems have been investigated and applied to control tannery effluents and mine discharges [5,6]. The problem of pollution following mine closure is particularly acute in many European countries where historical mining has been undertaken with poor pollution control and where the high costs of treatment make restoration a major issue [7].
This paper focuses on pollution in the Dhaleshwari River system in Dhaka, Bangladesh, where a new complex of tanneries at Savar was established in 2016 ( Figure 1). The Dhaleshwari river system is part of the Jamuna-Brahmaputra catchment and is the major source of freshwater to the Greater Dhaka Watershed. Flows in the river have decreased due to natural as well as human interventions, but augmentation of flow during the dry season in the Dhaleshwari and some other rivers is under consideration targeting pollution reduction and ecological restoration works planned under the Bangladesh River Masterplan and other initiatives. A new novel multi-branch model, INCA-Metals, has been developed to model pollution impacts in rivers and their tributaries, and the model has been applied to assess metal pollution and a range of restoration strategies. The INCA-Metals model builds on the existing suite of INCA models that have been developed over the past 20 years to simulate both terrestrial and aquatic systems [8][9][10][11]. The structure of INCA has been tested on many catchments, including 20 UK catchments, 25 catchments across Europe, and globally on large river systems such as the Mekong, the Ganges and the Volta [12][13][14]. Early major applications of INCA have been published in two special volumes of international journals [15,16]. The advantage of INCA is that the model seeks to incorporate the dominant mechanisms and processes operating so that a realistic and rapid assessment of environmental change can be evaluated. In this paper, a new version of the model, INCA-Metals, is described for rivers and catchments to simulate eight metals (cadmium, lead, zinc, mercury, arsenic, copper, chromium, manganese), as well as cyanide and ammonia. Cyanide and ammonia are included so that the metals model can be applied to mining waste issues. Cyanide is a key component of gold metal extraction processes, and ammonia is often released from the mining effluent treatment process. Both can be extremely toxic to river systems and hence have been included in the model. A new novel multi-branch model, INCA-Metals, has been developed to model pollution impacts in rivers and their tributaries, and the model has been applied to assess metal pollution and a range of restoration strategies. The INCA-Metals model builds on the existing suite of INCA models that have been developed over the past 20 years to simulate both terrestrial and aquatic systems [8][9][10][11]. The structure of INCA has been tested on many catchments, including 20 UK catchments, 25 catchments across Europe, and globally on large river systems such as the Mekong, the Ganges and the Volta [12][13][14]. Early major applications of INCA have been published in two special volumes of international journals [15,16]. The advantage of INCA is that the model seeks to incorporate the dominant mechanisms and processes operating so that a realistic and rapid assessment of environmental change can be evaluated. In this paper, a new version of the model, INCA-Metals, is described for rivers and catchments to simulate eight metals (cadmium, lead, zinc, mercury, arsenic, copper, chromium, manganese), as well as cyanide and ammonia. Cyanide and ammonia are included so that the metals model can be applied to mining waste issues. Cyanide is a key component of gold metal extraction processes, and ammonia is often released from the mining effluent treatment process. Both can be extremely toxic to river systems and hence have been included in the model.

The INCA-Metals Model
The Integrated Catchment Framework Model (INCA) is used as the main platform  for the new multibranch metals model. INCA is a catchment scale process-based model  to calculate pollutant transfer from diffuse sources and point sources to the catchment  outlet. To date, the INCA suite of models includes simulation of flow, nitrogen, phosphorus, sediments, chloride, carbon, pathogens, dissolved oxygen and mercury [8][9][10][11][17][18][19]. Many of the parameters and processes that have previously been successfully validated in other versions of INCA also apply to modelling of metals, such as the hydrology, landscape pathways, sediment settlement and entrainment, and flushing out components. Due to the dynamic nature of the model, variations in rainfall and temperature and changes of inputs, such as tannery treatment or mine discharges and runoff from diffuse pollution, can be investigated. This means the effects of climate change, population growth, land-use change and mitigation measures can be evaluated.
INCA-Metals is a mass-balance dynamic model that estimates the daily fluxes and concentrations of flow, ammonia, cyanide and eight metals in rivers. The eight metals are cadmium, lead, zinc, mercury, arsenic, copper, chromium and manganese. The flows and metals are simulated based on mass balance calculations accounting for contributions from various inputs and transformations. Key physico-chemical processes and stores are assumed in the landscape and within the river, and these are shown in Figure 1. The equations for INCA-Metals are based on those written for the Integrated Catchment model of Nitrogen (INCA-N [8-10]) but have been adapted to describe the metal adsorption to sediment, cyanide decay to ammonium and cyanide volatilization. The inputs and outputs are differentiated by landscape type and varied according to the environmental conditions, namely soil moisture and temperature. The model accounts for stocks of ammonium, cyanide and the eight metals in the soil and groundwater zones and in the stream reaches. The model also simulates the flow of water through the soil and groundwater zones from different soil and land use types to deliver mass to the river. This mass is then routed downstream after accounting for direct inputs from point sources, abstractions, within-river precipitation and sedimentation of metals, nitrification and cyanide decay and volatilization. Point sources such as discharges from tannery effluents or mine adits, waste dumps or tailings ponds, or effluent treatment facilities can be routed into any of the sub-catchments or reaches of the model.

The Simulation of Water Flow and Storage in the Landscape
The landscape mass balances of water, cyanide, ammonium and the eight metals are based on a 1 km 2 cell ( Figure 2). The inputs to the model can vary on a sub-catchment basis and according to soil or land-use type. These two factors allow the mass stored, process rates and hydrological pathways to vary spatially based on variations in soil moisture, temperature, adsorption potential and land management. The water volumes and the mass of cyanide, ammonium and the eight metals are summed based on the relative amounts of each land use or soil type within a sub-catchment and the output passed to the instream routing model. The hydrology component of the model is represented in Figure 3, where precipitation flows through the soil and groundwater zones, although there can be overland flow or direct runoff flow in high rainfall events.     The flow of water through the soil and groundwater zones is given by the following two equations: Soil Zone where q sz and q gz are the outflows from the soil and groundwater zones (m 3 s −1 km −2 ), p eff is the hydrologically effective rainfall (m 3 s −1 km −2 ), β is the base flow index (Ø) and T sz and T gz are the response times associated soil and groundwater zones (days). Within the soil zone, it is assumed the water can be partitioned into two volumes: drainage and retention. The drainage volume represents the water stored in the soil that responds rapidly to water inflow and drains under gravity; it may be thought of as macropore or drain flow: the flow that most strongly influences the rising hydrograph limb. The soil zone retention volume represents the water stored retained in the soil after gravity drainage; it responds more slowly than the drainage water and represents the majority of water in the soil.
The initial value of the soil zone drainage volume (V D , m 3 km −2 ) is calculated from a user-supplied initial soil zone flow (q sz,initial , m 3 km −2 ) and the soil water response time: where V D,initial is the initial soil water drainage volume (m 3 km −2 ), q sz,initial is the usersupplied initial soil water flow rate (m 3 s −1 ) and T sz is the user-supplied soil water time constant (days). The initial value of the soil zone retention volume (V R , m 3 km −2 ) is calculated from the supplied soil moisture deficit time-series (SMD, mm), an estimate of the maximum soil moisture deficit (SMD max , mm) and a parameter that describes the linear relationship between the soil moisture deficit and the soil zone retention volume, C 1 . The value of this parameter represents the ratio of the total water in the retention volume relative to the easily available water, namely that about the wilting point. The value of the parameter is derived through calibration and typically takes a value of 1 to 3. The value of SMD initial is estimated from the soil moisture deficit on day one of the simulation.
The groundwater initial volume (V gw , m 3 km −2 ) is estimated from the maximum size of the store and the proportion of pore space filled with water at the start of the model run: where V gw,initial is the initial groundwater volume (m 3 km −2 ), d eff,gw is the user-supplied maximum groundwater effective depth (m, active depth x effective porosity) and C 2 is the user-supplied proportion of filled pore-space (Ø). Given the complexity of the geology in most model applications, there is no attempt to separate the groundwater into specific yield and retention components; the former is the water able to drain from the rock under gravity, and the latter is the water retained against the pull of gravity. The soil zone drainage and retention volumes and the groundwater zone volume are recalculated at each time step to account for the changes in input and output.

The Simulation of the Transport, Storage and Transformations of Cyanide, Ammonium and Metals in the Landscape
The change in cyanide mass in the soil, m cn,sz (kg CN km −2 ) and groundwater, m cn,gz (kg CN km −2 ) stores are given by Equations (6) and (7): Soil Zone dm cn,sz dt = −m cn,sz q sz 86, 400 where C 3 , C 4 and C 5 are the rates of cyanide volatilisation and soil and groundwater zone conversion of cyanide to ammonium ( Figure 4). All other terms have been defined previously except S SMD , which is the soil moisture factor and describes the linear dependency of the rate of the soil zone processes on the soil moisture. It is assumed that there are no diffuse inputs of cyanide to the catchment. The first term on the right-hand side of Equation (6) represents the lateral transport of cyanide with the soil water to the stream, the second term represents cyanide volatilisation and the third term represents the decay of cyanide to ammonium. The first and second terms in Equation (7) represent the flux of cyanide into the groundwater from the soil zone and the lateral flow of cyanide with the groundwater into the stream; the third term represents the decay of cyanide to ammonium in the groundwater. The soil moisture factor is calculated at each time step as where SMD is the user-supplied daily soil moisture deficit time series (mm). The factor is scaled between 0 and 1 and describes the situation in which, as the soil dries, the rate of the soil processes declines. In addition, each process rate parameter is a function of the soil temperature where θ s is the soil temperature ( • C), C n is the soil process parameter and t q10 (Ø) and t Q10bas ( • C) are parameters determined by calibration. The parameter t q10 is the factor change in rate with a 10-degree change in temperature, and the parameter t Q10bas is the base temperature for the process at which the response is 1. The soil temperature is estimated from a seasonal relationship dependent on air temperature as follows where θ s is the air temperature ( • C) and C 6 is the maximum temperature difference between summer and winter ( • C). This relationship generates a seasonal pattern for each land use, which is controlled by the parameter C 6 . It is corrected for snow depth during winter months using Equation (4) from Rankinen et al. (2004) [20]. This temperature dependency is applied to all soil zone process parameters, not just those for cyanide. where C3, C4 and C5 are the rates of cyanide volatilisation and soil and groundwater zone conversion of cyanide to ammonium ( Figure 4). All other terms have been defined previously except SSMD, which is the soil moisture factor and describes the linear dependency of the rate of the soil zone processes on the soil moisture. It is assumed that there are no diffuse inputs of cyanide to the catchment. The first term on the right-hand side of Equation (6) represents the lateral transport of cyanide with the soil water to the stream, the second term represents cyanide volatilisation and the third term represents the decay of cyanide to ammonium. The first and second terms in Equation (7) represent the flux of cyanide into the groundwater from the soil zone and the lateral flow of cyanide with the groundwater into the stream; the third term represents the decay of cyanide to ammonium in the groundwater. The soil moisture factor is calculated at each time step as where SMD is the user-supplied daily soil moisture deficit time series (mm). The factor is scaled between 0 and 1 and describes the situation in which, as the soil dries, the rate of the soil processes declines. In addition, each process rate parameter is a function of the soil temperature = 10 − 10 10 where θs is the soil temperature (°C), Cn is the soil process parameter and tq10 (Ø ) and tQ10bas (°C) are parameters determined by calibration. The parameter tq10 is the factor change in rate with a 10-degree change in temperature, and the parameter tQ10bas is the base temperature for the process at which the response is 1. The soil temperature is estimated from a seasonal relationship dependent on air temperature as follows where θs is the air temperature (°C) and C6 is the maximum temperature difference between summer and winter (°C). This relationship generates a seasonal pattern for each land use, which is controlled by the parameter C6. It is corrected for snow depth during winter months using Equation (4) from Rankinen et al. (2004) [20]. This temperature dependency is applied to all soil zone process parameters, not just those for cyanide.

Ammonium-N
The change in ammonium mass in the soil, m nh4,sz (kg N km −2 ) and groundwater, and m nh4,gz (kg N km −2 ) stores are given by Equations (11) and (12): Soil Zone Groundwater Zone dm nh4,gz dt = βm nh4,gz q sz 86,400 V gw − m nh4,gz q gz 86,400 V gw − C 9 m nh4,gz 10 6 V gw (12) where C 7 , C 8 and C 9 are the rates of ammonium uptake by plants from the soil and the soil and groundwater ammonium nitrification rates ( Figure 4). The input mass of ammonium to the soil zone, m nh4,in (kg N ha day −1 ), represented by the first term on the right-hand side of Equation (11), includes livestock manure, fertiliser and wet and dry deposition. All other terms have been defined previously except S PGI , which is the plant growth index and describes the seasonal variations in the solar radiation and the subsequent control on plant growth during different seasons. The index is given by where C 10 is the day number associated with the start of the growing season. The second term on the right-hand side of Equation (11) represents the lateral movement of ammonium with flow from the soil zone to the stream, the third term represents the ammonium uptake by plants, the fourth term represents the conversion of ammonium to nitrate via nitrification and the fourth term is the increase in ammonium resultant from the decay of cyanide. The first and second terms on the right-hand side of Equation (12) represent the input of ammonium from the soil zone to groundwater and the lateral movement to the stream; the third term represents ammonium nitrification in the groundwater.

Metals
The following equations relate to each of the eight metals included in the model. Whilst the equations given here are generic to all the eight metals, within the model there are different parameters for each metal so that the transport, decay and storage of each are independent.
The change in mass of each metal in the soil, m metal,sz (kg km −2 ) and groundwater, m metal,gz (kg km −2 ) stores are given by Equations (14) and (15) Soil Zone Groundwater Zone dm metal,gz dt = βm metal,gz q sz 86,400 V gw − m metal,gz q gz 86,400 V gw − C 12 m metal,gz 10 6 V gw (15) where C 11 and C 12 are the rates of metal adsorption to the soil sediment and aquifer matrix. It is assumed that there are no diffuse inputs of metals to the catchment. The first term on the right-hand side of Equation (14) represents the lateral movement of metals transported with the flow to the stream. The second term represents the adsorption of a metal to the soil sediment. The first term on the right-hand side of Equation (15) represents the input of a metal from the soil to the groundwater zone by percolation, the second term represents the lateral flow of the metal from the groundwater to the stream and the third term represents the adsorption of metal to the aquifer matrix. The main flow paths for the metals are shown in Figure 5. represents the adsorption of metal to the aquifer matrix. The main flow paths for the metals are shown in Figure 5.

The Simulation of Water Flow and Storage in the River
The reach residence time constant, Treach (days) is calculated as where L is the reach length, qreach,out is the discharge from the reach and a and b are parameters relating the reach velocity to the discharge. The parameters a and b are determined through calibration, though b typically has a value of 0.67. The parameters can also be determined from tracer measurements. The change in the reach flow is calculated from the input-output mass balance of the form where qreach,in is the sum of the input flows from upstream, if not the headwater reach, point source effluent, diffuse inputs from the soil and groundwater zones and a loss via abstraction.

Cyanide, Ammonium, and Metal Process Equation: River System
In the river, the key processes are cyanide volatilisation and degradation, ammonium nitrification and metal loss due to sedimentation or precipitation. The reach mass balance includes the upstream water quality together with diffuse inputs from the soil and groundwater zones, as well as direct effluent discharges and abstractions ( Figure 6).
The cyanide mass, mcn,reach (kg) stored in a river reach is given by where the mass into the reach, mcn,reach_in (kg day −1 ) is the sum of the upstream input, point source effluent and the diffuse inputs from the soil and groundwater. The second term on the right-hand side of Equation (18) represents the mass transfer downstream with the flow of water, the third term represents the volatilisation of cyanide, the fourth term represents the decay of cyanide to ammonium and the fifth term represents any mass that is

The Simulation of Water Flow and Storage in the River
The reach residence time constant, T reach (days) is calculated as where L is the reach length, q reach,out is the discharge from the reach and a and b are parameters relating the reach velocity to the discharge. The parameters a and b are determined through calibration, though b typically has a value of 0.67. The parameters can also be determined from tracer measurements. The change in the reach flow is calculated from the input-output mass balance of the form dq reach,out dt = q reach,in − q reach,out T reach (17) where q reach,in is the sum of the input flows from upstream, if not the headwater reach, point source effluent, diffuse inputs from the soil and groundwater zones and a loss via abstraction.

Cyanide, Ammonium, and Metal Process Equation: River System
In the river, the key processes are cyanide volatilisation and degradation, ammonium nitrification and metal loss due to sedimentation or precipitation. The reach mass balance includes the upstream water quality together with diffuse inputs from the soil and groundwater zones, as well as direct effluent discharges and abstractions ( Figure 6).
The cyanide mass, m cn,reach (kg) stored in a river reach is given by where the mass into the reach, m cn,reach_in (kg day −1 ) is the sum of the upstream input, point source effluent and the diffuse inputs from the soil and groundwater. The second term on the right-hand side of Equation (18)  of water, the third term represents the volatilisation of cyanide, the fourth term represents the decay of cyanide to ammonium and the fifth term represents any mass that is removed from the reach by abstraction. The rate parameters describing cyanide volatilisation, C 13 , and decay, C 14 , are both temperatures dependent such that where θ w is the water temperature ( • C), which is assumed equal to the input air temperature.
The ammonium mass, m cn,reach (kg N) stored in a river reach is given by where the mass into the reach, m nh4,reach_in (kg day −1 ) is the sum of the upstream input, point source effluent and the diffuse inputs from the soil and groundwater. The second term on the right-hand side of Equation (20) represents the mass transfer downstream with the flow of water, the third term represents the nitrification of ammonium, the fourth term represents the increase in ammonium from the decay of cyanide and the fifth term represents any mass that is removed from the reach by abstraction. The rate parameters describing nitrification C 15 is temperature dependent.
The metal mass, m metal,reach (kg N) stored in a river reach is given by where the mass into the reach, m metal,reach_in (kg day −1 ) is the sum of the upstream input, point source effluent and the diffuse inputs from the soil and groundwater. The second term on the right-hand side of Equation (21) represents the mass transfer downstream with the flow of water; the third term represents the sedimentation of the metal.; and the fourth term represents any mass that is removed from the reach by abstraction. The rate parameters describing sedimentation C 16 is temperature-dependent.
where θw is the water temperature (°C), which is assumed equal to the input air temperature.
The ammonium mass, mcn,reach (kg N) stored in a river reach is given by where the mass into the reach, mnh4,reach_in (kg day −1 ) is the sum of the upstream input, point source effluent and the diffuse inputs from the soil and groundwater. The second term on the right-hand side of Equation (20) represents the mass transfer downstream with the flow of water, the third term represents the nitrification of ammonium, the fourth term represents the increase in ammonium from the decay of cyanide and the fifth term represents any mass that is removed from the reach by abstraction. The rate parameters describing nitrification C15 is temperature dependent. The metal mass, mmetal,reach (kg N) stored in a river reach is given by where the mass into the reach, mmetal,reach_in (kg day −1 ) is the sum of the upstream input, point source effluent and the diffuse inputs from the soil and groundwater. The second term on the right-hand side of Equation (21) represents the mass transfer downstream with the flow of water; the third term represents the sedimentation of the metal.; and the fourth term represents any mass that is removed from the reach by abstraction. The rate parameters describing sedimentation C16 is temperature-dependent. Within the soil zone, it is assumed that the water can be partitioned into two volumes: drainage and retention. The drainage volume represents the water stored in the soil that responds rapidly to water inflow and drains under gravity; it may be thought of as macropore or drain flow (i.e., the flow that most strongly influences the rising hydrograph Within the soil zone, it is assumed that the water can be partitioned into two volumes: drainage and retention. The drainage volume represents the water stored in the soil that responds rapidly to water inflow and drains under gravity; it may be thought of as macropore or drain flow (i.e., the flow that most strongly influences the rising hydrograph limb). The soil zone retention volume represents the water stored or retained in the soil after gravity drainage; it responds more slowly than the drainage water and represents the majority of water in the soil.
All the model constants and variables in the landscape equations and in the river equations are given in Tables 1 and 2.  Table 2. Constants and variables in the river mass-balance equations.

Symbol Definition Units
User-supplied inputs as constants

Tannery Pollution Globally and in Bangladesh
With about 300 million hides treated globally per year, leather tanning is one of the most environmentally problematic industries due to its high water consumption and toxic wastewater. Tannery effluents generally have a complex composition with very high levels of organic and inorganic pollutants as well as high concentrations of heavy metals, such as chromium, lead and cadmium. These can be extremely damaging to human health and aquatic ecosystems. Chromium tanning is the most common process due to its relatively short turnaround time and excellent quality yield. Around 90% of the world's leather gets Cr-tanned [21] using over 15,000 tons of Cr-salts annually. Only 20% of the chemical used in the production process is fixed, and therefore large fluxes of chromium are released via tannery wastewater into river systems. In Bangladesh, the tannery industry has recently been relocated from Hazaribag to Savar in Dhaka (Figure 7), with a total of 155 factories moved. This move was designed to relieve serious pollution in the Buriganga River and provide a central effluent treatment plant (CETP) at Savar. However, ongoing delays in effective operation of the CETP have called into question whether the pollution problem has been just transferred from the Buriganga River to the Dhaleshwari River ( Figure 7).

Impacts of Tannery Discharges on Pollution on Bangladesh Rivers
Tanneries generate significant polluting discharges into river systems unless effective effluent treatment is applied. Table 3 shows the mean chemistry from 16 published studies of tannery waste in Bangladesh. There are high levels of pollutants in the effluents, and Table 4 shows the average of these, with very low dissolved oxygen levels, corresponding with high BOD (Biological Oxygen Demand) and COD (Chemical Oxygen Demand) concentrations as well as dissolved solids discharges. Nutrients such as nitrate and phosphorus are also significant, as is chloride (or salt content). Most concerning are the metal concentrations, with extremely high lead, zinc, cadmium and chromium concentrations.
The high concentrations of metals, nutrients and organic compounds found in the tannery effluents have been shown to be deleterious to the aquatic eco-system. There are a wealth of studies on the harmful impacts of the tannery discharges on the Buriganga River water quality, sediments and fish populations [22]. These hazards can also harm human health through bio-magnification with metals such as chromium causing acute renal failure, gastrointestinal haemorrhage, lung cancer, asthma and nervous system damage, as shown in Table 5. Other metals such as lead, for example, can impair children's development, cause memory loss and insomnia, cause vomiting, diarrhoea and kidney damage (Table 5).

Impacts of Tannery Discharges on Pollution on Bangladesh Rivers
Tanneries generate significant polluting discharges into river systems unless effective effluent treatment is applied. Table 3 shows the mean chemistry from 16 published studies of tannery waste in Bangladesh. There are high levels of pollutants in the effluents, and Table 4 shows the average of these, with very low dissolved oxygen levels, corresponding with high BOD (Biological Oxygen Demand) and COD (Chemical Oxygen Demand) concentrations as well as dissolved solids discharges. Nutrients such as nitrate and phosphorus are also significant, as is chloride (or salt content). Most concerning are the metal concentrations, with extremely high lead, zinc, cadmium and chromium concentrations.
The high concentrations of metals, nutrients and organic compounds found in the tannery effluents have been shown to be deleterious to the aquatic eco-system. There are a wealth of studies on the harmful impacts of the tannery discharges on the Buriganga River water quality, sediments and fish populations [22]. These hazards can also harm human health through bio-magnification with metals such as chromium causing acute renal failure, gastrointestinal haemorrhage, lung cancer, asthma and nervous system damage, as shown in Table 5. Other metals such as lead, for example, can impair children's development, cause memory loss and insomnia, cause vomiting, diarrhoea and kidney damage (Table 5).

An Experiment on the Dhaleshwari River System and the Savar Tannery Complex
The pollution in Hazaribag has been of national concern since the 1980s, when the government first instructed tanneries to control their pollution footprint [40,41]. In 2001, The High Court of Bangladesh ordered the tanneries to apply the necessary measures to limit their toxic discharges within one year, following the submission of a written petition by the Bangladesh Environmental Lawyers Association (BELA). In 2002, the government approved a relocation plan to move the tanneries to Savar district (adjacent to the Dhaleshwari river), and in 2003 the Bangladesh Government signed an industry transfer agreement with the Small and Cottage Industries Corporation (BSCIC) and major leather tannery industry associations. Due to a series of political and economic issues, the relocation only formally took place in 2017. While currently over 90% of the tanneries have now moved to Savar, the new effluent treatment plant is still incomplete and not operating fully and efficiently. The introduction of tannery discharges in Savar has created additional pollution in the Dhaleshwari River system. Studies have revealed a gradual deterioration of its water quality over the past two decades [42][43][44][45], with several urban and industrial discharges along the river. Heavy metals have risen in water, sediments, and with alarming levels recorded in fish, which have been designated "unfit for human consumption" [43].
Few studies have been carried out to monitor and model pollution scenarios in the Dhaleshwari. Two studies assessed the relocation of the tanneries to Savar, but only focused on the impact on the Buriganga. Islam (2018) [40] used the Water Analysis Simulation Program (WASP) model to investigate the impact of the tanneries' relocation on the physicochemical parameters of the Buriganga. BOD, ammonia, nitrate and phosphate levels decreased significantly, but DO levels did not improve. Whitehead et al. [41] studied heavy metals (Cd, Cr, As, Cu) and treatment scenarios in the Buriganga River using the INCAmetals model. The simulations revealed that unless an 85% wastewater treatment level was attained, the heavy metals would still be above permissible limits in the river system.
In order to assess the pollution impacts of the new Savar tannery complex, a sampling and chemical analysis programme was conducted in August and September 2019. Figure 8 shows the sampling sites that are located upstream and downstream of the Savar tannery discharge point. The sampling strategy involved sampling at a depth of 0.6 m from water surface, using a Teflon-coated plastic bottle (autoclavable and acid-washed) with samples transfer in iceboxes during fieldwork, then in a refrigerator (at 4 • C) before their shipment to Oxford, UK. Chemical analysis was conducted using Inductively Coupled Plasma-Mass Spectrometry (ICP-MS). Measurements by ICP-MS were undertaken using a PerkinElmer NexION 350D ICP-MS at the Department of Earth Sciences, University of Oxford. Each metal was calibrated from a series of calibration standards, which were automatically produced by the coupled Elemental Scientific prepFAST M5 syringe-driven autosampler.
Samples were taken weekly for six weeks from all the sites (Figure 8), and chemical analysis conducted. Table 6 shows the average concentrations of metals in the upstream and downstream sites and the percent change between the two sites. The data show a consistent increase in pollution from upstream to downstream with chromium and cadmium, showing 254% and 65% increases, respectively. The metals concentrations at the upstream site are significant in themselves and this represents pollution from upstream discharges. So the new Savar tannary discharges are just exacerbating the existing pollution in the river system. The concentrations in the river are quite low because the river flows were high during the sampling programme and hence river dilution was significant. It should be noted that the sampling period coincided with the peak production season (post-Eid Al Adha), which follows a lunar calendar. As EID occurs at different times of the year it is important to consider the pollution levels in different seasons and especially during low flow periods when dilution is low. Therefore, in order to assess the impact of the river under low flow conditions, the river's water quality was further investigated under low-flow conditions. Based on the historical flow dataset of the Dhaleshwari River system, Mimouni [39] has estimated the dilution factor to be of the order of 130. This implies that metal concentrations in the river would be approximately 130 times greater in the dry season compared to the wet season (Table 6). Of course, other environmental variables change on a seasonal basis such as temperature, and it is important to bear these in mind during any modelling study. As described above, the INCA equations are temperaturedependant and so will change with changing temperature. The model also considers the changing water velocities and hence residence times so that decay and sedimentation rates will vary from high-flow to low-flow conditions.  Samples were taken weekly for six weeks from all the sites (Figure 8), and chemical analysis conducted. Table 6 shows the average concentrations of metals in the upstream and downstream sites and the percent change between the two sites. The data show a consistent increase in pollution from upstream to downstream with chromium and cadmium, showing 254% and 65% increases, respectively. The metals concentrations at the upstream site are significant in themselves and this represents pollution from upstream discharges. So the new Savar tannary discharges are just exacerbating the existing pollution in the river system. The concentrations in the river are quite low because the river flows were high during the sampling programme and hence river dilution was significant. It should be noted that the sampling period coincided with the peak production season (post-Eid Al Adha), which follows a lunar calendar. As EID occurs at different times of the year it is important to consider the pollution levels in different seasons and especially during low flow periods when dilution is low. Therefore, in order to assess the impact of the river under low flow conditions, the river's water quality was further investigated under low-flow conditions. Based on the historical flow dataset of the Dhaleshwari River system, Mimouni [39] has estimated the dilution factor to be of the order of 130. This implies that metal concentrations in the river would be approximately 130 times greater in the dry season compared to the wet season (Table 6). Of course, other environmental variables change on a seasonal basis such as temperature, and it is important to bear these in mind during any modelling study. As described above, the INCA equations are temperature-dependant and so will change with changing temperature. The model also considers the changing water velocities and hence residence times so that decay and sedimentation rates will vary from high-flow to low-flow conditions.  Since other polluting sources also release their effluents into the Dhaleshwari, attributing the pollution loads to the tanneries' discharges remains challenging. In order to identify the heavy metals that were particularly impacted by the tanneries' effluent during the sampling period, an amalgamation of different factors/conditions (spatial and temporal) that corroborate the potential contribution of the CETP effluent to the heavy metals concentrations recorded were taken into consideration. Any impact attributed to the CETP effluent was also based on their effect on the Dhaleshwari's water quality vis-à-vis international/regional surface water quality standards (Table 7) and the link of each analysed heavy metal to the tanning process. Strong impacts of three heavy metals were recorded: Cr, Cd and As. Considering the environmental threats that high levels of these heavy metals entail, the improvement of CETP's treatment performance is imperative to ensure the compliance of its Cr, Cd and As discharges. This must be prioritized during all seasons not to exhaust the river's assimilative capacity and to minimize the toxicity threats to the wider ecosystem and local population.

Modelling the Savar Discharge and Mitigation Strategies
In order to determine the percentage level of upgrade needed in the CETP and highlight the impact of its improved performance on the water quality of the Dhaleshwari, different improved CETP treatment scenarios were simulated using INCA-Metals. This is the first modelling study to have targeted the impact of the relocated tanneries' discharges at Savar on heavy metals in the Dhaleshwari River. This has involved setting up the process-based dynamic flow and quality model INCA Metals for the Dhaleshwari River catchment so that we could simulate a series of river reaches from its source down to the Buriganga. As shown in Figure 8, we model reaches of the river that correspond to the sampling locations and then add additional reaches down to the confluence point with the Buriganga River (Figure 7). Knowing the flow data for the river and the effluent discharge flow rate, then, by mass balance, we can estimate the effluent discharge concentrations. Using these concentrations, we can then evaluate the model performance, and the model simulates the concentrations of a range of metals with some accuracy, as shown in Figure 9. We can then use the model to simulate the impact of the discharge on the river system under a range of seasonal conditions and also a range of mitigation scenarios. Figure 10 shows a set of simulations along the river under low flow (i.e., dry season/worse case) conditions together with a set of treatment scenarios for the three heavy metals identified to have strongly impacted the Dhaleshwari during the sampling period. The concentrations along the river were simulated under the estimated low flow conditions to further highlight the gravity of the situation. It can be seen that levels of Cr, Cd and As are very high (black line) and well above international environmental standards if no treatment improvement is made (Table 7). These three key metals are of particular interest because of their impact on human health and fisheries. The results suggest that the CETP at Savar needs to operate at much higher levels of efficiency to control the river pollution, depicted by the other lines on the graphs representing improvements to the treatment plant under a range of scenarios from 65% to 95% reductions in discharge concentrations. In order to meet reasonable environmental conditions, an improvement of approximately 85% is required. Of course, a big issue with water quality modelling is that of uncertainty as decay parameters have not been measured explicitly but assumed from other studies [7,41]. Thus, the modelling results are subject to uncertainty and may vary along the river system under different assumptions. Another limitation is the lack of information regarding other discharges along the river, which were consequently not incorporated into the model used. Thus, the treatment improvement percentages suggested in our study are non-conservative and may vary, depending on many assumptions, but are nevertheless important in order to highlight the level of pollution control needed to safeguard the water quality of the Dhaleshwari.
is made (Table 7). These three key metals are of particular interest because of their impact on human health and fisheries. The results suggest that the CETP at Savar needs to operate at much higher levels of efficiency to control the river pollution, depicted by the other lines on the graphs representing improvements to the treatment plant under a range of scenarios from 65% to 95% reductions in discharge concentrations. In order to meet reasonable environmental conditions, an improvement of approximately 85% is required. Of course, a big issue with water quality modelling is that of uncertainty as decay parameters have not been measured explicitly but assumed from other studies [7,41]. Thus, the modelling results are subject to uncertainty and may vary along the river system under different assumptions. Another limitation is the lack of information regarding other discharges along the river, which were consequently not incorporated into the model used. Thus, the treatment improvement percentages suggested in our study are non-conservative and may vary, depending on many assumptions, but are nevertheless important in order to highlight the level of pollution control needed to safeguard the water quality of the Dhaleshwari.

Policy Recommendations
There are several clear messages for policy recommendations resulting from the modelling and monitoring analysis in this paper. The Savar plant may require a treatment upgrade from the current process solely based on a biological treatment to one that incorporates a chemical primary treatment to better handle the large concentrations of heavy metals released in the tanneries' wastewater [39]. Such an upgrade could lead to significant improvements in water effluent quality. Operational management must be ensured

Policy Recommendations
There are several clear messages for policy recommendations resulting from the modelling and monitoring analysis in this paper. The Savar plant may require a treatment upgrade from the current process solely based on a biological treatment to one that incorporates a chemical primary treatment to better handle the large concentrations of heavy metals released in the tanneries' wastewater [39]. Such an upgrade could lead to significant improvements in water effluent quality. Operational management must be ensured to achieve a much higher treatment efficiency, and further attention should be reserved for the tanneries' sludge management, which also impacts both the land and water quality of the area. New technology should be considered that can treat the particularly difficult waste from tanneries [46]. The Leather Working Group, a private multi-stakeholder group responsible for auditing and certifying leather producers, should consider adjusting their requirement, for tanneries to ensure wastewater treatment of their effluents, to the political, social and economic context of Bangladesh. It might be proactive to incorporate water-quality-oriented auditing of the Savar treatment plant to prompt action for improved control. For the government, it will be necessary to establish initiatives to enhance the uptake of clean production practices. The Department of Environment (DoE) as the main regulatory agency should liaise with other stakeholders to establish a realistic, time-bound approach to achieve SDG 6.3. Additionally, the DoE could establish ambient water quality standards based on the different pollutants released by local industries and the wastewater treatment capacity in Dhaka. It would be advisable to upgrade the capacity of the DoE to include real-time monitoring of effluent discharges and river water quality, as well as continuously ensuring public dissemination of water quality data. Long-term goals should target the analysis of riverine sediment and biota. The INCA-Metals model covers the whole of the Dhaka region, and this could form the basis of a real-time system to forecast pollution events and identify hotspots for decision-makers.

Conclusions
This paper provides the first modelling analysis of water quality changes in the Dhaleshwari since the relocation of tanneries from Hazaribag to Savar. It introduces the INCA-Metals model and its novel application in the Dhaka watershed. The paper then introduces the background of the study. It explains our field findings of increases in heavy metals concentrations from the upstream to downstream sites from the tanneries' discharge point, expected to reach serious levels during the upcoming dry season, particularly for Cr, Cd and As. The results suggest that an 85% reduction in discharge concentrations is required to meet internationally approved standards for these particular metals. It should be emphasized that these findings need further investigation during peak and off-peak periods throughout the year in order to validate the results. The failure of the CETP to achieve sufficient pollution reduction calls into question the future viability of relying on central technical interventions, rather than considering other forms of recycling, reuse, or reduction of heavy metals in the process of leather tanning. For Bangladesh to achieve sufficient progress toward SDG 6, particularly target 6.3, addressing heavy metals pollution from leather tanneries will be a crucial step.
Effluent treatment facilities are considered a key step in ensuring reduced pollution of river systems from individual factories or from a cluster or sites, as is the case at Savar. At the moment, it appears that the treatment plant at Savar is underperforming. Sustainable Development Goal (SDG) 6 [1] for water and sanitation addresses the key issues of water quality and quantity. As in many countries, Bangladesh is striving to progress towards SDG 6.3 to ensure the adequate supply of good quality water for many uses including public supply, industrial supply and to restore rivers to good water quality and ecological status. Moving whole industrial sectors to new locations without putting in place working treatment plants is counterproductive as it merely pollutes yet another stretch of the river system. Urgent action needs to be taken to improve the treatment facilities at Savar. The findings in this paper offer useful insights for other studies evaluating metals transport or the impact of particular industries upon the aquatic environment. Future work may usefully consider additional issues including speciation of metals, seasonality and other variables impacting water quality change over time and mechanisms for implementing improved monitoring systems.