Numerical Study of Sulfonamide Occurrence and Transport at the Near-Shore Area of Laizhou Bay

: Antibiotics are extensively applied in aquaculture for the treatment of microbial infections and to improve production. Among various antibiotics, sulfonamides (SA) are popular in ﬁsh farming, and SA residues in the aquatic environment have detrimental e ﬀ ects on both ecosystem and human health. Understanding the fate of SA in the aquatic environment is a basic necessity to provide an approach for solving the current problem. In this study, a two-dimensional lattice Boltzmann model was introduced to investigate the transport and occurrence of SA in Laizhou Bay, a prosperous aquaculture area in China. The model is based on the shallow-water equations and advection-di ﬀ usion equation with the Bhatnagaar–Gross–Krook scheme. Experimental data are used to verify the model after numerical simulations and the results illustrate the accuracy of the proposed model. This model provides a potential universal method for the simulation of the fate of antibiotics in the aquatic environment.


Introduction
Aquaculture is the farming of aquatic organisms by individuals or groups with the goal of enhancing their production. Rapid population growth in the 21st century requires a sharp increase in seafood production [1]. Seafood provides a good source of high-quality protein and forms an important cash crop in many parts of the world. However, the majority of aquaculture production relies on the use of antibiotics and input of other artificial substances to treat or prevent disease and to accelerate growth. Therefore, studies on antibiotic residues in the aquatic environment as well as in fish products are numerous [2][3][4]. These residues could increase the occurrence of antibiotic-resistant bacteria in the aquatic environment, increase antibiotic resistance in fish pathogens, lead to variations of the bacterial flora both in sediments and water systems, and transfer these resistance determinants to human pathogens and bacteria of terrestrial animals [5].
Approximately 70 classes of antibiotics have been reported in the aquatic environment [6]. Among these, sulfonamides (SA) are the most widely used in aquaculture cultivation [2,4]. Furthermore, SA have been detected in the muscles of fish species through bioaccumulation according to recent research [7]. Given sufficient time, SA can accumulate in the body and adversely impact human health. Moreover, bacteria will become resistant to SA, which tends to contribute to multidrug and cross resistance. It is difficult to predict future negative effects caused by the potential environmental and health risks. Therefore, the fate of SA in the aquatic environment, as a basis of the current problem, plays an important role in relevant research fields.
With regard to the investigation of antibiotic occurrence and transport, many scholars focus on empirical methods. The most popular method is liquid chromatography tandem mass spectrometry (LC-MS/MS) [8][9][10]. However, most empirical methods focus on the occurrence, not the transport process, of antibiotics. With respect to the numerical method, there are many mathematical models for pollutant transport in the aquatic environment [11,12], while very few attempts to simulate the transport of antibiotics have been made [13]. Therefore, a promising numerical method, the lattice Boltzmann method (LBM), is used in this study to address this issue. The LBM was developed during the last few decades. It evolved directly from the Lattice Gas Automata (LGA) introduced by Frisch [14]. The LBM is recommended as a mesoscopic method due to its nature. It can describe particle microscopic behaviors with a simpler form and accurately reflect flow movements at a macroscopic level compared with traditional methods. Consequently, it has been extensively used by researchers for the simulation of fluid flows [15][16][17]. This study is the first to study the transport of antibiotics using the LBM.
In this paper, the shallow-water equations and depth-averaged advection-diffusion (AD) equation are solved using the LBM, and then are used to simulate the transport and occurrence of SA in Laizhou Bay, a highly polluted sea area of China. The results of this study can be used to aid the understanding the environmental risks of antibiotic mixture exposures to aquatic species and human health in the future.

Shallow Water Equations
The two-dimensional (2D) shallow water equations (SWE) including the continuity equation and momentum equation [18], are: where the subscripts i and j are the space direction indices and the Einstein summation convention is used; x j is the Cartesian coordinate, taking x, y and z in turn; u j is the velocity component which takes u and v corresponding to that in x and y and directions, respectively. h is the water depth; t is the time; ν is the kinematic viscosity; Z b is the bed height of the datum plane and F i is the force term and can be regarded as bed shear stress τ bi , which is calculated by: where C b is the bed friction coefficient, which is calculated by C z = h 1 6 /n b , where C z is the Chezy coefficient given by the Manning coefficient n b .

Depth-Averaged Advection-Diffusion Equation
In general, the 2D-AD equation is: where c is the solute concentration per unit water depth; K d is the degradation coefficient; and K ij is the diffusion coefficient.

LBM for Shallow Water Equations
Based on the LBM with a D2Q9 lattice (Figure 1), the scheme with the enhanced force term [15,16] can be expressed as: where is the weight factor: = 4/9 for α = 0; = 1/9 for α = 1,3,5,7; = 1/36 for α = 2,4,6,8; is the distribution function of particles; is the local equilibrium distribution function; Δt is the time step; τ is the single relaxation time; is the velocity vector of a particle in the α link. The external force term can be written as: For the D2Q9 lattice, each particle moves one lattice at its velocity only along the eight links indicated by 1-8; "0" indicates resting particles with zero speed. The particle velocity is: where e = Δx/Δt; Δx is the lattice size. The local equilibrium distribution function is determined by: Furthermore, the water depth h and velocities can be calculated by following equations: The viscosity is defined as: The external force term F α can be written as: For the D2Q9 lattice, each particle moves one lattice at its velocity only along the eight links indicated by 1-8; "0" indicates resting particles with zero speed. The particle velocity is: where e = ∆x/∆t; ∆x is the lattice size. The local equilibrium distribution function f eq α is determined by: Furthermore, the water depth h and velocities u i can be calculated by following equations: The viscosity is defined as:

LBM for Depth-Averaged Advection-Diffusion Equation
The LBM equation for the AD equation can be expressed as: g eq α is the local equilibrium distribution function; τ s is the single relaxation time for advection-diffusion equation; and s α is the velocity vector of a particle in the α link which can be described as: where s = ∆x/∆t. The local equilibrium distribution function g eq α is determined by: in which λ ij is the non-dimensional coefficient defined by: Therefore, after calculating the h and u i via the LBM from the hydrodynamic model, the concentration c can be calculated by:

Numerical Study
Based on the work of [8], two types of SA, including sulfamethoxazole (SMX) and sulfamethazine (SMZ), were simulated with the coupled model. Then, a comparison between simulation results and the measured data of near-shore locations P1, P2, and P3 ( Figure 2) was made to verify the proposed model.

Research Area
Laizhou Bay (Figure 2) is one of the three largest bays of the Bohai Sea, north of China. The bay is an important fishery base in Shandong province, providing a substantial volume of aquatic products. Currently, with the development of urbanization and industrialization, aquaculture in Laizhou Bay is prosperous. Moreover, more than 20 rivers flow into the bay. Therefore, antibiotic residues from industries are inevitably flushed into the bay through riverine inputs.

Acquisition of the Data
GIS was used to obtain the topographic data of Laizhou Bay. The accurate tidal data in this study were calibrated by the previous research [19]. The experimental comparison data were derived from [8].

Boundary Conditions
The coastlines were treated with wet-dry boundary condition [20], whereas the open boundary was set as a tidal wave [19].

Parameter Setting and Sensitivity Analysis
The research area was divided into 140 × 230 grids with ∆x = ∆y = 460 m in consideration of the required accuracy and computational efficiency of the proposed model. ∆t was 4 s according to the ∆x/∆t ≫ , which satisfies the stability condition. We used τ = 0.6 for the LBM-SWE equation and τ = 1.01 for the LBM-AD equation. The viscosity was calculated by Equation (10). The friction coefficient Cb is in the range of 0.61-22.7, which depends on the different hydrodynamic conditions. The half-life of SA is = 18.48 ℎ [21], and the degradation process can be expressed as: The initial still-water depth was set to 20 m. The initial antibiotic concentrations of local rivers are shown in Table 1, and these values were placed at the grids of the nine estuaries as pollution sources.

Acquisition of the Data
GIS was used to obtain the topographic data of Laizhou Bay. The accurate tidal data in this study were calibrated by the previous research [19]. The experimental comparison data were derived from [8].

Boundary Conditions
The coastlines were treated with wet-dry boundary condition [20], whereas the open boundary was set as a tidal wave [19].

Parameter Setting and Sensitivity Analysis
The research area was divided into 140 × 230 grids with ∆x = ∆y = 460 m in consideration of the required accuracy and computational efficiency of the proposed model. ∆t was 4 s according to the ∆x/∆t u, which satisfies the stability condition. We used τ = 0.6 for the LBM-SWE equation and τ s = 1.01 for the LBM-AD equation. The viscosity was calculated by Equation (10). The friction coefficient C b is in the range of 0.61-22.7, which depends on the different hydrodynamic conditions. The half-life of SA is t 1 2 = 18.48 h [21], and the degradation process can be expressed as: The initial still-water depth was set to 20 m. The initial antibiotic concentrations of local rivers are shown in Table 1, and these values were placed at the grids of the nine estuaries as pollution sources. Diffusion coefficients (Kij) varied with different solutes and external environmental factors, ranging from 0.02 to 28 m 2 /s. Therefore, it was necessary to perform the sensitivity test on Kij in consideration of model stability and fitting of experimental data (Figure 3). Isotropic diffusion is applied in this study, i.e., Kxx = Kyy and Kxx = Kyy = 0. When Kxx = Kyy = 10 m 2 /s, the simulation results were basically consistent with the experimental data [8]. When Kxx, Kyy were reduced by 10% and 20%, the simulation results decreased evenly 9.97% and 20%, respectively. Likewise, when Kxx, Kyy were raised by 10% and 20%, it also leads to an average increase of 10.03% and 20% in the simulation, respectively.  Diffusion coefficients (Kij) varied with different solutes and external environmental factors, ranging from 0.02 to 28 m 2 /s. Therefore, it was necessary to perform the sensitivity test on Kij in consideration of model stability and fitting of experimental data (Figure 3). Isotropic diffusion is applied in this study, i.e., Kxx = Kyy and Kxx = Kyy = 0. When Kxx = Kyy = 10 m 2 /s, the simulation results were basically consistent with the experimental data [8]. When Kxx, Kyy were reduced by 10% and 20%, the simulation results decreased evenly 9.97% and 20%, respectively. Likewise, when Kxx, Kyy were raised by 10% and 20%, it also leads to an average increase of 10.03% and 20% in the simulation, respectively.  Figure 4 shows that the SMX and SMZ concentrations change over time at the three selected points. SMX is the dominant SA compared to SMZ and poses a high risk to aquatic organisms. Out of all the rivers, Xiaoqing River and Guangli River are the two most contaminated; therefore, the SA concentrations at P1 and P2 are higher than that at P3. The transport process also indicates increasing pollution of the aquatic environment without control of antibiotic abuse.  Figure 4 shows that the SMX and SMZ concentrations change over time at the three selected points. SMX is the dominant SA compared to SMZ and poses a high risk to aquatic organisms. Out of all the rivers, Xiaoqing River and Guangli River are the two most contaminated; therefore, the SA concentrations at P1 and P2 are higher than that at P3. The transport process also indicates increasing pollution of the aquatic environment without control of antibiotic abuse.

Results and Discussion
Comparison between simulation results and experimental data [8] with different SA and locations after 42 h are shown in Table 2. The results illustrate that the proposed model is sufficiently accurate to obtain the occurrence of SA. The results also emphasize that the concentrations of SA are high in the near-shore area-especially for SMX at P2, which is comparatively high in the marine area at 80 ng/L. Comparison between simulation results and experimental data [8] with different SA and locations after 42 h are shown in Table 2. The results illustrate that the proposed model is sufficiently accurate to obtain the occurrence of SA. The results also emphasize that the concentrations of SA are high in the near-shore area-especially for SMX at P2, which is comparatively high in the marine area at 80 ng/L.  Table 2, the simulation value of SMZ at P3 is available with 0.087 ng/L compared to the undetected experimental measurement. Owing to limitations of instrument precision, the results are deficient at small magnitudes. Therefore, the numerical method, the lattice Boltzmann method is recommended for this work owing to its simplicity and accuracy. Less than 5 h are required to simulate 42 h of real transport, highlighting the efficiency of the model. In this case, the flows have Re < 250, belonging to the laminar flow, a range within which the LBM has very good stability. Furthermore, complex boundary conditions can be processed well, and parallel operations can be conducted.
The occurrence of SMX and SMZ has been plotted in Figure 5 to show the spatial distribution of SA after 30 days. On the last day, the pollution covers the whole bay, and SA concentrations > 0.005 g/L are dominant.  In Table 2, the simulation value of SMZ at P3 is available with 0.087 ng/L compared to the undetected experimental measurement. Owing to limitations of instrument precision, the results are deficient at small magnitudes. Therefore, the numerical method, the lattice Boltzmann method is recommended for this work owing to its simplicity and accuracy. Less than 5 h are required to simulate 42 h of real transport, highlighting the efficiency of the model. In this case, the flows have Re < 250, belonging to the laminar flow, a range within which the LBM has very good stability. Furthermore, complex boundary conditions can be processed well, and parallel operations can be conducted.
The occurrence of SMX and SMZ has been plotted in Figure 5 to show the spatial distribution of SA after 30 days. On the last day, the pollution covers the whole bay, and SA concentrations > 0.005 g/L are dominant. Due to the high density of human activities and the many rivers flowing into the bay, the high SA concentrations in the model show the dominant impact of riverine inputs on the marine environment. This finding coincides with previous studies of the sources of antibiotics in marine systems [8,22]. Most estuaries are located in the northwest of Laizhou Bay, leading to higher pollution levels in this area. The bay is semi-enclosed with poor water exchange, resulting in high-level Due to the high density of human activities and the many rivers flowing into the bay, the high SA concentrations in the model show the dominant impact of riverine inputs on the marine environment. This finding coincides with previous studies of the sources of antibiotics in marine systems [8,22]. Most estuaries are located in the northwest of Laizhou Bay, leading to higher pollution levels in this area. The bay is semi-enclosed with poor water exchange, resulting in high-level environmental risk. Moreover, the results for the selected points are lower than those in most river estuaries, which may result from seawater dilution where rivers discharge into the bay [23].

Conclusions
This paper proposes a two-dimensional model, coupling hydrodynamic and water quality models. The LBM is used to discretize the model and the bed shear stress is also considered in this study. Moreover, the coefficients for the diffusion and relaxation times are calibrated, and the coupled model is applied to simulate the transport and occurrence of SA, including SMX and SMZ in Laizhou Bay. The numerical results are in good agreement with the experimental data, showing that the proposed model is able to predict the transport of sulfonamides accurately. It is reasonable to expect that the proposed model, as an efficient tool, has the potential to be used to investigate antibiotic transport problems in complex environmental conditions.