Integration of numerical models to simulate 2D hydrodynamic/water quality model of contaminant concentration in Shatt Al-Arab River with WRDB calibration tools

Abstract: The hydrodynamic model is essential for building a water quality model for rivers, lakes, estuaries, and other water systems. Most model software, such as HEC-RAS, can perform a complex hydrodynamic surface water body and limitations to represent water quality for the corresponding area. In contrast, other models, like WASP, can simulate a wide range of contaminants in a multidimensional geometry of rivers, estuaries, lakes, and reservoirs. Still, it requires flow information from separate hydrodynamic models. This article aims to develop a comprehensive water quality model of the Shatt Al Arab River south of Iraq by linking HEC-RAS with WASP. A variety of software techniques has sequentially been used. This software includes GIS for DEM modification, HEC-RAS for the hydrodynamic model, Python code with PyCharm to run the external coupler, WASP software for advective and dispersive contaminant transport, and finally, WRDB software for full calibration process and results display. The results showed successful transportation of flow information had been achieved. Moreover, the article described an effective calibration process by plotting comparison graphs and statistical summaries to make the appropriate decision. Another goal of this work is to collect the equations and associated reaction rates of source/sink kinetic for eutrophication’s state variables.


Introduction
The modeling of surface water is complex and developing [1]. Professionals disagree with the "optimal" way of modeling rivers, lakes, estuaries, and coastal waters [2]. Even after more than a century, the basic approach to surface water modeling has not altered because all models are built on three fundamental principles: mass, momentum, and energy conservation [3,4]. The most general classification of modules used in surface waters modeling is analytical or numerical models [1]: An analytical model has an exact mathematical solution to the governing equations that describe processes in a body of water. An analytical solution example for estimating the concentrations of dissolved oxygen (DO) along rivers is the Streeter-Phelps (1925) equation [5]. A numerical model is a discretized representation of a mathematical equation system that explains processes in a water body. The computing domain is discretized into cells, and the partial differential governing equations are approximated by a set of algebraic equations solved by the iteration or the matrix inversion method. In recent years, computer simulation techniques have gained popularity in scientific studies, particularly regarding studies of the aquatic environment. Many computer models have been created and are effectively used in many countries today [1]. Review research studies from the Web of Science focused on the fate and transport of water quality modules in waterbodies, such as EFDC, CE-QUAL-W2, WASP, Delft3D, AQUATOX, and MIKE, are listed in Table 1. The highest models utilized in uses and citations were EFDC and CE-QUAL-W2. The United States and China were the most frequent users of such models [6]. Fu et al. [7] searched the Scopus database for 50,530 publications on water quality models published between 1935 and 2018. 76% of the publications (38,542) were published between 2003 and 2018. Most of these articles were published by authors with locations in the United States, followed by China, the United Kingdom, Germany, Canada, and Australia.
The following paragraphs will cover a variety of research that used standard computer models with their most prominent results.
Ramos-Ramírez et al. [8] undertook water quality simulations of the dispersion of Chrome III (Cr III) in the Bogotá River, Cundinamarca, Colombia. The WASP model was employed in the simulation's construction to reflect the Cr III dispersion under the research circumstances. Using the most recent 1D, 2D, and 3D modeling software, Marvin and Wilson [9] described the thorough hydrodynamic modeling done for the estuary of Salmon River and its branches as a portion of an overall flood hazard analysis. Alam et al. [10] investigated the impacts of rising temperatures and solar radiation from climate change on the Sitalakhya River's water quality using WASP. Chuersuwan et al. [11], with the help of WASP, aimed to empower local authorities on water quality management in the Lamtakhong river's basin, Thailand. Yen et al. [12] investigated pollution sources and water quality comprehensively, followed by implementing the WASP model to develop water quality requirements in Carp Lake, Taiwan. In surface water quality modeling, the model coupling is a significant trend. The integrated modeling technique is more sophisticated and difficult to master than the conventional. The main obstacles to the coupling process are a lack of observed data, bathymetry information, model parameter estimates, and computational time [13]. WASP has been proved in the literature to be capable of simulating river transport. It has been integrated with several hydraulic tools, frequently for water quality modeling via external Linkage to hydraulic models for water flow information. The prominent instances are as follows.
Chueh et al. [14] simulated the copper concentration using SWAT and WASP in the Erren River. SWAT was used to compute non-point source (NPS) soil erosion, while WASP was used to simulate copper concentrations in both the aqueous and sediment phases. Zhu et al. [15] investigated a unidirectional link of the urban hydrological and hydraulic models CS-TARP, WASP, and EFDC that were used to simulate the influence of sewer overflows on nutrients in the Chicago River. Define et al. [16] combined the Regional Ocean Modeling System (ROMS) with WASP to conduct a thorough examination of water quality in Barnegat Bay, New Jersey. A WASP model integrated with the HEC-HMS model estimated mercury transport and fractionation in the Xiaxi River, China [17]. Yin and Seo [18] developed a three-dimensional model using an integrated hydrodynamic Environmental Fluid Dynamics Code (EFDC) and water quality model WASP for the Ara Canal in Korea that connects the Han River and the Yellow Sea. SWAT and WASP simulations of the watershed NPS load and internal flow in the Cedar Creek Reservoir, according to Ernst and Owens [19], had the greatest effect on chl a and TP. Wool et al. [20] employed WASP and hydrodynamics data from EFDC to estimate eutrophication, nutrient circulation, and DO kinetic in the Neuse River estuary, United States. Rodriguez and Peene [21] developed a two-dimensional hydrodynamic, transport, and water quality model in order to determine the allowable total maximum daily load (TMDL) for oxygendemanding materials. The 2D vertically averaged hydrodynamic model (EFDC) was combined with WASP.
This article aims to develop a comprehensive water quality model of the Shatt Al-Arab River south of Iraq by linking HEC-RAS with WASP. Moreover, the article described an effective calibration process by plotting comparison graphs and statistical summaries to make the appropriate decision. Another goal of this work is to collect the equations and associated reaction rates of source/sink kinetic for eutrophication's state variables.

Description of the study area
The Tigris and Euphrates Rivers meet near Al-the Qurna district in southern Iraq to form the Shatt Al-Arab River [22] (Figure 1). The Shatt Al-Arab River is a 192-km-length tidal river that flows south-eastwards, passing through Basrah and then discharging into the Arabian Gulf [23]. The River's width changes along its length, from 250 to 300 m near the Euphrates-Tigris confluence to 600 m around the town of Basrah and 2,000 m at the estuary [24]. For the last 95 km of its course, the River forms part of the border between Iraq and Iran [25]. In addition to transportation, the Shatt Al-Arab River has an important role in delivering water for domestic use, irrigation, and manufacturing [26]. Several water treatment plants along the River divert water for household use. Currently, contributions of River feeders have been decreased due to surrounding governments' policies, which resulted in a significant increase in TDS values in the Shatt Al-Arab River due to the effect of the Arabian Gulf [27]. Data for the year 2014 were used to develop the model. These data include cross-sections, daily flow and hourly stage hydrographs, temperature, conductivity, and TDS. Table 2 shows the location of the sampling stations along the River and the contribution of each data in the model building. A reported discharge from a separate model and a recorded hourly water stage for 2014 extended from 1 February up to 30 June were chosen as input upstream and downstream BC, respectively.

Basic model elements
The movement of the River has an impact on the physical transport and transmission of contaminants. Water motions at various scales and forms substantially impact the aggregation and distribution of nutrients, dissolved gases, temperature, and the distribution of microorganisms and plankton. As a result, the water quality model's accuracy is highly dependent on the hydrodynamic model's description of circulation. Therefore, three submodels were used sequentially to integrate the Shatt Al Arab River modeling: the hydrodynamic model of HEC-RAS, the heat model of WASP, and the eutrophication model of WASP.

Hydrodynamic model (HEC-RAS)
The study area, Shatt Al-Arab River, can be considered a shallow water surface because the width-to-depth ratio is more than 10 [28,29]. Therefore, the depth-averaged 2D hydrodynamic model is adopted to simulate the hydrological parameter variation in plain view (longitudinal and transverse direction) in HEC-RAS. However, computing the 2D Saint Venant equations [30,31] demands greater processing power, resulting in longer run durations. Furthermore, numerical instability of the equations can occur in areas of the 2D mesh where the flow path or depth profile changes abruptly. A refinement and shorter step will be required to avoid an unstable model [32].

Water temperature model (WASP)
WASP, developed by EPA, can simulate a wide range of contaminants and temperatures in a multidimensional geometry of rivers, estuaries, lakes, and reservoirs. One of the essential physical features of surface waterways is water temperature [33]. It is important in hydrodynamic and water quality research [2]. The temperature may be thought of as a measure of how much heat energy is kept in a water volume, and equation (1) describes the heat energy and temperature relationship [34]: where H is the heat content, V is the volume, ρ is the density of water, and C p is the specific heat of water   [33] or equilibrium temperature and surface heat exchange coefficients. Water temperature sources and sinks related to surface heat exchange can be expressed as follows: where A s is the surface area in m 2 , t is the time in the day, and H n is the net thermal energy flux (Watt/m 2 ) where H s is the shortwave solar radiation flux, H a is the longwave radiation flux, H e is the evaporation heat loss flux, H c is the heat conduction flux, H sr is the reflected short wave solar radiation, H ar is the reflected long-wave radiation, and H br is the back radiation from the water surface, also see Figure 2 for more explanation. Stephan Boltzmann's fourth-power radiation law is utilized for H br [29,34] s a (5) Shortwave solar radiation penetrates the surface and decays exponentially with depth according to Bears law: To compute the heat flux components, the related parameters with typical values are listed in Table 3. Indeed, some models compute surface fluxes and regard bed energy transfer as insignificant [33].
The equilibrium temperature technique assumes that the net surface heat exchange is zero. This definition of the equilibrium temperature allows describing H n as follows: where K aw is the coefficient of surface heat exchange, T w is the water surface temperature, and T e is the equilibrium temperature.

Eutrophication model
The core principle behind constructing a mass balance equation for a waterbody is to consider all materials entering and exiting the waterbody through direct material addition, such as runoff and loading, advection and  dispersion transport processes, and physical, chemical, and biological changes [35]. The 2D mass transport equation is

HC t HuC x
HvC y where H is the water depth, C is the vertical average concentration of the water quality constituent, t is the time, u and v are longitudinal and lateral advective velocities, respectively, E x and E y are longitudinal and lateral diffusion coefficients, respectively, S L is the direct loading rate, S B is the boundary loading rate, S K is the kinetic transformation rate; the source is positive, and the sink is negative. Depending on the types of state variables representation, different degrees of complexity have been recognized from simple Streeter-Phelps to nonlinear DO balance. Eight state variables that can contribute to eutrophication are included in Table 4 with their notations. It should be noted that all abbreviations appearing in the following equations have been explained in Appendix.

Phosphorus
Inorganic phosphorus (IP) is consumed by phytoplankton for growth and assimilated into the biomass of phytoplankton. Endogenous respiration and nonpredatory mortality return phosphorus from the phytoplankton biomass as organic phosphorus (OP) in either dissolved or particulate form. OP is mineralized into IP at a dependency temperature rate [36].

Nitrogen
The nitrogen composition kinetics is similar to that of the phosphorus cycle. NH 3 -N and nitrate are consumed by phytoplankton and assimilated into phytoplankton biomass for growth. The rate at which each is absorbed is controlled by its concentration concerning the total amount of inorganic nitrogen (NH 3 -N plus nitrate) available. Through endogenous respiration and nonpredatory death, nitrogen is recovered from the phytoplankton biomass pool as dissolved and particulate organic nitrogen ON and NH 3 -N. ON is mineralized to NH 3 -N at a temperature-dependent rate, and NH 3 -N is then transformed to nitrate at a temperature-and oxygen-dependent nitrification rate. In the absence of oxygen, nitrate may be transformed into nitrogen gas at a temperature and oxygendependent denitrification rate.

Carbonaceous biochemical oxygen demand (CBOD)
In addition to man-made sources and natural runoff, detrital phytoplankton carbons generated by algal death are major sources of CBOD. The principal CBOD loss mechanisms are oxidation, settling, and denitrification. When phytoplankton dies, organic carbon is produced that can be oxidized. Recirculation of phytoplankton carbon to CBOD is represented by first-order expression and stoichiometric oxygen-to-carbon ratio 32/12.

Phytoplankton
The kinetics of phytoplankton plays a major role in eutrophication, impacting all other state variables. Where phytoplankton increases due to photosynthesis and are lost via respiration, death, and settling, Figure 3 ( where S k4 is the reaction term, P is the phytoplankton population, G P1 is the growth rate constant, D P1 is the death plus respiration rate constant, and K S4 is the settling rate constant. Figure 4 depicts an overview of all transformation processes of the EUTRO system with related state variables.

Models coupling and development
WASP does not support multidimensional hydrodynamics simulation, but its adaptable structure allows importing data from an external hydrodynamics model. External hydrodynamic models coupled to WASP through linkage files can anticipate transport information for complicated water bodies. Before the transformation process begins, it should be verified for quality. Consequently, this check is discussed in detail in the findings section.

External linkage
The transport information output of the hydrodynamic model, such as volumes, flows, velocity, etc., is reported in the result file, which was then extracted and utilized by WASP [37]. To make this connection properly, an external coupler written in Python v.3.7 ( Figure 5), developed by [38], is used to link the two models via the development of an Application Program Interface (API). The various information in space and time that the hydrodynamic model sends to WASP contains segments' volume, depth, velocity, temperature, and salinity.

Segment description
The segment provides specific geometry information. To properly transmit data between the two model grids, the WASP segments were reproduced to HEC-RAS 2D cell boundaries by the hydrodynamic linkage file. It is preferable in modeling the mesh size of the HEC-RAS model (implicit finite-difference model) must be finer than the WASP model (explicit finite-difference model) to achieve convergence during the analysis [16]. In another ward, to prevent numerical errors for WASP, the HEC-RAS 2D computational grid cells should be collected for use in WASP ( Figure 6).

Model calibration
Tracers, whether synthetic or natural, such as salinity, chloride, dye, or heat, are frequently employed to calibrate the dispersion coefficients for a model network. A chemical tracer is a nonreactive (conservative tracer) chemical conveyed through a body of water. The first step in modeling more complicated water quality factors is to set up and calibrate a tracer. Five common statistical error methods  have been adopted in Water Resource Data Base (WRDB) calibration to give an adequate view for accepting or refusing the model outputs. These methods were coefficient of determination (r 2 ), mean absolute error (MAE), root-mean-square error (RMSE), normalized RMSE (N RMSE), and index of agreement (d) [39][40][41][42][43].
3 Results overview and discussion

HEC-RAS results
Calibration of HEC-RAS 2D output for water surface elevation is done by changing Manning's roughness coefficient (n). Three distinct values of "n" were considered: n = 0.02, n = 0.025, and n = 0.035; then chose, an appropriate "n" value that gives allowable criteria ranges between simulated and observed water stages. As it appears from Table 5, there is an acceptable convergence between the calculated and measured values at n value is equal to 0.025.  Figure 8. This figure illustrates the depth variation affected by tidal phenomena at different locations along the Shatt Al Arab River. These graphs were constructed by WRDB software. Now, WASP is ready to be initiated for contaminants simulation.

EUTRO model results
A 2D concentration distribution is displayed on the window of the WRDBGIS software as a colored spatial grid through   Figure 11 shows TDS concentration variation through a tidal period at the downstream end of the River, where the value increased gradually during the rising phase to extend the high concentration value for approximately 10 km from the downstream end, then decreased through the falling phase. After proceeding with the modeling process, it is seen that the penetration distance is increasing with time, where at date and time (2/16/2014 16:00), it is found to be 22.4 km with a TDS range of 9,000-10,000 mg/L (Figure 12a), and at date and time (3/17/2014 03:00), it is found to be 28.7 km with a concentration range of 8,000-9,000 mg/L (Figure 12b      Another way to display the WASP results is an x/y plot as a time series variation for the constituents. It is possible to choose several segments (cells) along the river to show the variation of TDS concentration at different stations, as shown in Figure 13. From this figure, the concentration of TDS increased toward the estuary. At the Qurnah station, the TDS ranged from 500 to 1,500 and from 1,500 to 2,000 mg/L at the Basrah station, while at the Sehan station was 1,500-3,000 mg/L. Moreover, one can illustrate the transversal distribution by selecting segments that lie on

Conclusion
This article aims to develop a comprehensive water quality model of the Shatt Al Arab River south of Iraq by linking HEC-RAS with WASP. A variety of software techniques has sequentially been used. This software includes GIS for DEM modification, HEC-RAS for the hydrodynamic model, Python code with PyCharm to run the external coupler, WASP software for advective and dispersive contaminant transport, and finally, WRDB/WRDBGIS software for full calibration process and results display. Therefore, three submodels were used sequentially to integrate the Shatt Al Arab River modeling: the hydrodynamic model of HEC-RAS, the heat model of WASP, and the eutrophication model of WASP. Data for the year 2014 were used to develop the model. These data include cross-sections, daily flow and hourly stage hydrographs, temperature, conductivity, and TDS. A reported discharge from a separate model and a recorded hourly water stage for 2014 extended from 1 February up to 30 June were chosen as input upstream and downstream BC, respectively. The study area, Shatt Al-Arab River, can be considered a shallow water surface because the width-to-depth ratio is more than 10; therefore, the depth-averaged 2D hydrodynamic model is adopted. Eight state variables (NH 3 , NO 3 , PO 4 , PHYT, CBOD, DO, ON, and OP) that can contribute to eutrophication are included with their equations of kinetic reaction. The results indicate a good flow path transformation through the River network. The calibration tool can give a clear vision to make the right decision about the model's accuracy.

Conflict of interest:
The authors state no conflict of interest.
Data availability statement: Most datasets generated and analyzed in this study are included in this submitted manuscript. The other datasets are available on reasonable request from the corresponding author with the attached information.   Preference for ammonia uptake -