Simulation of pollutant diffusion in vegetation open channel based on LBM-CA method

This paper proposes a pollution diffusion model that accurately assesses changes in instantaneous river pollution in vegetation open channels. The model is established based on cellular automata and lattice Boltzmann method (LBM-CA). Flow influence coefficients are incorporated into cellular automata (CA) to represent the effect of vegetation on pollutant diffusion, while the lattice Boltzmann method (LBM) is utilized to simulate flow in vegetation open channels and obtain the flow influence coefficients for each cellular. The results show that the LBM-CA model has high accuracy and that pollutants tend to accumulate in vegetation areas, thereby extending the residence time of pollutants. The model incorporates pollution limits, allowing the prediction of basin pollution levels at specific times. The LBM-CA model provides a method for simulating pollutant diffusion in natural rivers.


Introduction
In natural rivers, predicting pollutant diffusion is important. Convection diffusion equations are usually used to solve pollutants' diffusion behavior in nonvegetated rivers. However, there is a lot of submerged vegetation in natural rivers, and the existence of vegetation changes the distribution of river flow velocity (Jia et al. 2022). The convection diffusion equation is difficult to solve. Therefore, it is necessary to establish a model that can simultaneously simulate the velocity of a vegetated open channel and the prediction of pollutant diffusion. This is to achieve accurate monitoring and prediction of river pollution degrees.
As mentioned earlier, most water quality models are based on the convection diffusion equation. Models can be categorized into 1-D, 2-D, and 3-D models based on their dimensions. Streeter and Phelps (1925) proposed the first water quality model, a simple 1-D oxygen balance model. Thompson et al. (2004) used the MIKE 11 model to simulate wetland hydrology. Hadgu et al. (2014) used QUAL2K model to simulate nonpoint source pollution in the catchment area of the Darugu River, and the results show that the model can better reflect the measured data. The above models are all 1-D models. Since the 1-D model can only simulate the longitudinal water quality diffusion, it cannot  (1992). Debele et al. (2008) integrated SWAT and CE-QUAL-W2 water quality models, and the results showed that the model can be used to assess and manage water resources in complex basins composed of highland basins and downstream water bodies. In recent years, 3D river simulation is also developing rapidly. Jin et al. (2007) combined EFDC model, Reynolds average N-S equation, and SAV model to propose LOEM 3-D lake water quality model. Karamouz and Mahani (2021) simulated the discharge of the coastal sewage treatment plant through a Delft-3D model to assess the impact of flood on water quality. The above models usually need to solve partial differential equations simultaneously, which leads to complex calculations and low efficiency. It is difficult to simulate water quality rapidly. Therefore, researchers explored other deterministic methods of water quality simulation. Among them, cellular automata (CA) is a new model in recent years. Simple local rules make CA show great plasticity, which makes it possible to use it to build model systems in many cases. Karafyllidis (1997) first used the important tool of GIS spatial analysis "cellular automata" to simulate the diffusion and drift of water pollutants. However, its model structure is too simple, and its basic data are too abstract. This results in distortion of the spatial position, shape, and concentration field distribution information of the pollution diffusion in the analysis results. Wang et al. (2009) established a spatiotemporal model of water pollution zone movement and diffusion based on the principle of CA. Chukwuma (2016) evaluated the irrigation water quality index of rivers through CA. Abhijith and Mohan (2020) constructed a 1-D chlorine reaction transport model based on CA according to random walking particles. Bai et al. (2020) studied the response mechanism of nonpoint source pollutant load to land use change and different precipitation scenarios using CA and HSPF models. However, no researcher has used CA to simulate pollutant diffusion in vegetation open channels. The evolution rules of CA do not consider the influence of vegetation.
In the simulation of pollutant diffusion, the velocity of the river basin is an important parameter. In the vegetation open channel, due to the existence of vegetation, the velocity distribution of the flow will change greatly. In order to better embed velocity information into the CA, the lattice Boltzmann method (LBM) is selected to simulate the vegetation open channel. The advantages of LBM, such as simplicity, efficiency, and easy treatment of boundary conditions, in simulating fluid flows have been demonstrated. LBM is a mesoscopic model with advantages in describing flow details. LBM can more realistically describe vegetation characteristics. Moreover, the LBM has relatively short simulation time and low complexity. LBM is also established on the basis of CA. The two methods have inherent advantages in subsequent integration. So, this paper chooses LBM to simulate the ecological open channel. Chapman and Cowing (1939) proved that in the case of low Mach number, Euler equation describing ideal fluid and Navier-Stokes equation describing viscous fluid can be regarded as zero-order and first-order approximations of Boltzmann equation, respectively. Since then, the LBM has developed rapidly in the field of computational fluid dynamics. Currently, there are several researchers who have applied the LBM to the simulation of vegetation open channels. Gac (2015) used the large eddy simulation method based on LBM algorithm to calculate the vertical velocity profile of open channel flow with submerged and nonsubmerged vegetation. Yang et al. (2017)

Cellular automata
Cellular automata (CA) is a lattice dynamics model with discrete time, space and state, and localized spatial interactions and temporal causality. It simulates the space-time evolution of complex systems. The state of the cellular changes according to certain rules, and its evolution rules are local. The cellular automata affect the changes of surrounding neighbors' step by step through local changes. In the end, it spreads to every part of the system through time accumulation and space iteration.

Evolution rules of cellular automata in vegetation open channel
To simulate pollutant movement and diffusion in the vegetation open channel, the flow body area can be divided into the same square unit. Assuming that each square cell of the water surface is a cellular, each cell size is ∆ (Fig. 1). It is found in the pretest that pollutants' mass decreases gradually from the center to the periphery. The initial distribution of pollutants basically conforms to the normal distribution. Therefore, the normal distribution is adopted as the initial state of the instantaneous point source. Moreover, the diffusion form of pollutants in the flow body is close to a circle, so the Moore-type neighbor cellular model is used for simulation in this paper (Griffeath and Moore 2003).
The diffusion of pollutants can be expressed by Fick's law (Fick 1855) where J is the diffusion flux, D is the diffusion coefficient, C is the volume concentration of the diffusion substance, and C x is the concentration gradient. This law proposes that the diffusion of substances always spreads from high concentration to low concentration.
According to Fick's law, the diffusion movement of pollutants in water is mainly affected by many factors such as flow velocity, flow direction, and diffusion speed. The constructed cellular space ( Fig. 1) needs to include the above influencing factors. So the state of the cellular C t i,j at time t can be expressed as where i is the row number and j is the column number; Water indicates whether the cellular is located on the side wall or the flow body. If the cellular is a sidewall region (solid), Water = 0 adopts the evolution rule of the sidewall. If the cellular is flow body (liquid), then Water = 1 adopts the evolution rule of flow body. M t i,j is the pollutant mass of the cellular at time t. w t i,j is the wind speed. u t i,j is the flow velocity. wd t i,j and ud t i,j is the direction of wind and flow. Figure 2 shows the weight of pollutants in the central cellular and neighbor cellular at time t. First, consider the diffusion of the central cellular and the surrounding north, south, west, and east axes. As a first-order approximation,  the weight of diffused pollutants should be proportional to the weight difference of pollutants between cellular (Karafyllidis 1997). Therefore, diffusion in these four directions can be expressed by the following equation: where m is the still diffusion coefficient. The diffusion amount of pollutants between diagonal cellular should be less than that between neighbors' cellular. The weight of pollutants diffused from the central cellular to the diagonal cellular can be expressed by the following equation (Wang et al. 2009): where d is the oblique diffusion coefficient. The sum of outflow and inflow of pollutants in the cellular cannot be greater than the total amount of original pollutants in the cellular, so m + md ≤ 0.25 . The most accurate simulation effect can be obtained when m = 0.084 and d = 0.16 (Karafyllidis 1997).
Velocity, wind, and direction were not considered in the previously established evolution rules. That is, Eq. (4) is the diffusion rule of pollutants in still water simulation (Fig. 3). In moving water, in addition to the diffusion simulation of pollutants, it is also necessary to simulate the movement of pollutants. Pollutant movement simulations can be divided into movement with flow and movement with wind. This experiment is conducted indoors, so the wind influence is ignored. According to the research results of Burn and Mcbean (1985), in the flow body, the evolution rule of pollutant diffusion needs to add the flow influence coefficient W on the basis of Eq. (4). The equation is as follows: where u t i,j is the flow velocity of the central cellular at time t and u max is the maximum flow velocity in the flow.
Therefore, the diffusion evolution rules of pollutants in vegetation open channels can be expressed as Due to flow velocity, in addition to the diffusion of pollutants, it is also necessary to consider the state of pollutants drifting in water. Since diffusion and water transport are carried out simultaneously, the time required for each diffusion step is equal to one-time step, and the time step for each diffusion is where t CA is the time step, Δ CA is the size of the cell, and v P is the diffusion velocity of pollutants. The further drift distance S can be calculated according to the following equation: It can be found in Eqs. (6) and (8) that the flow velocity at each point will have a certain impact on pollutant diffusion. For pollutant diffusion in vegetation channels, (a) Initial state (b) One step diffusion fluid mechanics researchers have put forward many famous formulas, such as the longitudinal dispersion coefficient model proposed by Fischer (1967), but the velocity distribution needs to be input into this model. On the basis of Fisher's research, Seo and Baek (2004) used beta density formula to describe the lateral distribution of flow velocity in the river and then calculated the longitudinal dispersion coefficient. Although this method can successfully predict the transverse distribution law of the flow velocity and obtain the longitudinal dispersion coefficient by three times integration, the relevant parameters in the beta density formula need to use some actual measured data, including the value and position of the maximum flow velocity. Shiono and Knight (1991) proposed that the lateral distribution of river flow velocity can be obtained by solving the N-S equation integrated along the water depth. Wang and Huai (2016) simplified the solution process of the longitudinal dispersion coefficient by using this method and combining Fourier series changes. However, in summary, it can be found that the basic pollutant diffusion models need the flow velocity as the input parameter. Therefore, velocity information needs to be obtained before CA simulation. In this paper, the LBM is used to simulate the flow of the vegetation open channel. The specific content will be described in the "Lattice Boltzmann method" section.

Boundary conditions of cellular automata in vegetation open channel
In the numerical model, another important problem is the model boundary conditions. This paper assumes that the pollutants in the flow will not be adsorbed after encountering the side wall of the open channel and will all be bounced back to the flow. Therefore, the channel boundary adopts a fixed boundary to meet the preset requirements. Taking the north boundary as an example, the equations are as follows:

Lattice Boltzmann method
Lattice Boltzmann method originated from cellular automata, and CA and LBM have inherent advantages in later integration. Therefore, LBM is used to simulate the vegetation open channel in this paper. The theory holds that fluid velocity can be dispersed into many particles. The total number of particles with a certain velocity in a space element at a certain time and position is called the particle distribution function. Every particle in space migrates in a given direction and collides and integrates at the lattice nodes. In LBM, if particles collide with each other, the total number of particles in this space will also change. The rate of change of particle distribution function from initial state to final state Ω is called the collision operator. The average time interval between two collisions is called relaxation time ( ). During the collision movement of particles, the particle distribution function gradually approaches the equilibrium distribution function ( f eq ), which is called relaxation in physics.

Governing equation of the lattice Boltzmann method in vegetation open channel
This paper uses the D2Q9 model for modeling. Figure 4 is a particle velocity direction distribution diagram of the D2Q9 model. The discrete speed and weight are shown in Table 1.
In the LBM, the motion of the particle is described by the Boltzmann equation where f i x i t is the particle velocity distribution function along direction i at time t moment x, ∆ is the length of a lattice, and t is the time step. In this paper, we use the classical LBGK collision model, whose collision function is expressed as (Bhatnagar et al. 1954) (10) As mentioned before, during the collision movement of molecules, the value of the particle distribution function will gradually approach the equilibrium distribution function. The expression of the function is (Chen and Doolen 2003) The LBM can be split into two steps: collision and streaming. The collision step is as follows: The streaming steps are as follows: The relationship between fluid viscosity and relaxation factor is For the LBM, mass and momentum conservation is also required in the simulation. The sum of distribution functions on each lattice is the macroscopic velocity of the fluid The momentum can be seen as the weighted value of the average velocity c i on the lattice and the distribution function The parameter selection in simulation is discussed in the "Parameter selection in ecological open channel LBM simulation" section. (14)

Boundary conditions of the lattice Boltzmann method in vegetation open channel
In LBM, it is an important problem to accurately simulate boundary conditions. Corresponding boundary conditions are given for the 2-D open channel in the east, west, north, and south directions.
• West boundary: nonequilibrium rebound approach In practical applications, the velocity components at the boundary are usually known, such as the inlet velocity of the channel flow. In this paper, we adopt a method proposed by Zou and He (Zou and He 1996) to calculate three unknown distribution functions perpendicular to the boundary under equilibrium conditions. In order to reduce the simulation error, the inlet lattice flow velocity for this simulation is set to 0.1. • East boundary: fully developed approach The length of the open channel can ensure the full flow of the fluid, so the full development scheme is adopted for the east boundary. This method considers that when the flow in the channel reaches full development, the physical quantities such as density and velocity will not change in the mainstream direction. The corresponding boundary processing format is implemented as follows: • North boundary/south boundary: bounce-back approach Since the north-south boundary is a solid boundary (no slip boundary), the standard rebound scheme is adopted for both boundaries. The rebound scheme refers to that the incident particles at the solid boundary will bounce back to the reflow volume domain (Fig. 5). In this paper, the bounce-back approach is also applied to the obstacle surface. The particle velocity on the surface and inside the obstacle is 0, and the specific equation will not be repeated. When programming, it should be noted that there are still four corners at the west boundary and the east boundary that are not constrained. These four corners need to be individually constrained in combination with the conditions of each boundary; otherwise, divergence may occur.

Model validation
In order to verify the correctness of the LBM model, the test data of Han (2014) are used in this paper. The channel used in the test is a rectangular open channel. The channel's length is 11.0 m, width 0.3 m, and height 0.45 m. The measuring point is located 6 m downstream of the water intake. Figure 6 shows the difference between measured and simulated values. In Fig. 6, the X-axis is y∕b , and the Y-axis is u∕u * . Both coordinates are dimensionless values. The frictional flow velocity can be expressed by the following equation: where g is gravity acceleration, H is channel water depth, and S is channel slope. The black lines in Fig. 6 represent simulated values, and the scattered points represent measured values. It can be clearly seen in Fig. 6 that the scattered points are evenly scattered around the simulation line. Figure 7 shows the error between measured and simulated values. The X-axis is the experimental value, and the Y-axis is the simulation value. In the middle of the open channel, the difference between the measured value ( u e ) and the simulated value ( u s ) is very small, with an error of within 5%. The error of all measuring points is within 10%, with high consistency. The above results demonstrate the reliability of the established LBM model.

Experiment
In this paper, three trials of open channel flow are used to verify the model. The three trials (Table 2) are the channel with scattered vegetation, the channel with aggregated vegetation (Zong and Nepf 2010), and the compound channel with vegetation.
The dispersed vegetation channel experiment was carried out in the hydraulic public test hall of China Agricultural University. The channel's length is 6.3 m, width 0.8 m, and height 0.6 m. Simulated aquatic plants were used to simulate vegetation in the experiment. The height of single vegetation is 0.06 m and the diameter is 0.03 m. The vegetation is arranged equidistant, and the vegetation spacing is 0.06 m × 0.06 m (Fig. 8a).
The vegetated compound channel was carried out in China Institute of Water Resources and Hydropower Research. The channel's length is 10.0 m, width 1.0 m, and height 1.0 m. The left side of the channel is a smooth channel, and the right side is covered with vegetation. The experiment adopts the real vegetation Cymbidium tenuifolia. The height of a single plant is 0.06 m, the diameter is 0.03 m, and the vegetation zone is 0.010 m × 0.075 m evenly arranged (Fig. 8b). Put rhodamine tracer at the Sect. 0.5 m away from the inlet of the channel. Contaminants are injected instantaneously under lateral pressure. Measurement points are set up upstream and downstream for the test. Rhodamine concentration was measured using YSI (rhodamine probes, YSI Inc., Yellow Springs, OH). The first YSI is placed at x = 2.5 , and the second YSI is placed at x = 5.5 . Set test points (yellow points in Fig. 9) at z = 0 and z = 0.25 , respectively, for two working conditions. The 3-D Ultrasonic Doppler (ADV) current meter was used to measure the velocity. ADV can measure the main flow velocity and vertical velocity of the flow field in real time. ADV is installed on the guide rail that can be moved horizontally and vertically above the channel. In order to minimize the impact of ADV noise, the sample size of each measuring point in this experiment is about 2000 times, and the experimental time is 200 s. To filter the overall data, the test uses 5 times the sample standard deviation as the interval.

Dispersed vegetation open channel simulation
The test point is located in the middle of the vegetation area. In this paper, flow velocity data under 80 m 3 /h and 100 m 3 /h flows are selected to verify the LBM model. All data in Fig. 10 are dimensionless. Figure 10 is a comparison diagram of LBM simulated value and measured value. It can be seen from Fig. 10 that the flow velocity is large near both sides of the channel and small in the middle vegetation area. The longitudinal flow velocity fluctuates periodically along the transverse direction. The greater the flow velocity, the greater the velocity fluctuation. Figure 11 plots the error between simulated and measured values. Most errors between simulated and measured values are within 10%. This indicates that the LBM can effectively simulate flow in dispersed vegetation open channels. Some of the points in the figure deviate relatively far, possibly due to the presence of vegetation leading to the formation of vortex streets in the water flow. The flow rate fluctuates periodically. Figure 12 shows the flow and turbulent kinetic energy (TKE) distribution clouds of the open channel at a flow of 100m 3 /h. Due to vegetation, the distribution of main flow has changed, and the main flow area tends to the side without vegetation. Behind each vegetation, there are separate dead flow zones of a certain length (Han et al. 2014). Also, it can be seen that the higher the flow velocity, the higher the TKE. The maximum

Aggregated vegetation open channel simulation
In order to verify the accuracy of LBM model in simulating ecological open channel on longitudinal velocity, the test data of Zong and Nepf (2010) are used in this paper. The diameter D of the aggregated vegetation community is 0.42 m and the diameter d of the cylinder is 0.006 m (Fig. 13). The vegetation is arranged staggered, and the density of aggregated vegetation can be described by solid volume fraction ( = Nd 2 ∕D 2 ). In this simulation, the inlet condition uniform flow U 0 is 0.098 m/s. The model established by LBM is consistent with the experiment. The calculation domain is a rectangular channel with a length of 13 m and a width of 1.2 m. The rigid vegetation community simulated by cylindrical group is arranged on the longitudinal central axis of the channel, and its upstream front end is 3 m from the entrance. In the experiment, the velocity of middle water depth is approximately replaced by the average velocity along the water depth. Therefore, this simulation is a 2-D plane of middle water depth. It does not consider the 3-D effects caused by the free water surface and boundary conditions of the channel bottom. The simulated flow and TKE are shown in Fig. 14. It can be seen from Fig. 14 that flow diversion begins about 1D away from the vegetation community. When flows pass through vegetation, the velocity on both sides increases. When x = 2/D, the velocity on both sides of the vegetation community reaches the maximum. The TKE is also large, which indicates strong energy exchange. Vegetation community has small velocity and small TKE, which indicates that vegetation has an obvious flow resistance effect (Xiang et al. 2019). Unlike the dispersed vegetation channel, the dead flow zone occurs behind the vegetation community and the length of the dead flow zone is longer. This reflects the influence of two different vegetation forms on flow. After the dead water zone, the flow velocity returns to steady flow. Figure 15 shows the comparison curves between the simulated results of longitudinal velocity u (along y = 0) and transverse velocity v (along y = D/2) and the experimental values. It can be seen from the diagram that the transverse flow velocity begins to decrease near the vegetation community at y = D/2 and then rises rapidly and reaches its maximum value at 0 < x/D < 1. The change trend of transverse (a) Δ=0.02m =0.51 (b) Δ=0.005m =0.80  Fig. 14b, which means that the flow diversion is the strongest at the maximum transverse flow velocity. Then, along the x direction, the transverse flow velocity v decreases rapidly and reaches a steady state gradually.
To estimate the LBM reliability, this paper also uses Fluent to conduct relevant simulations. In Fluent simulation, the geometric conditions of the model are the same as those of LBM. The model is meshed using unstructured grids. In this paper, the standard k − turbulence model is used for simulation. The water flow inlet adopts a velocity inlet boundary, and the water flow outlet adopts a pressure outlet boundary. Vegetation and sidewalls are defined as nonslip boundary conditions. Figure 15 shows the simulation results.
It can be seen in Fig. 15 that the flow velocity trend of LBM simulation value is basically in agreement with the experiment value. There are two main differences between simulated and experimental values. One is that near upstream of the vegetation community, the transverse flow velocity tends to decrease. This part of the change is not shown in the experiment, but both LBM and Fluent simulations confirm that the transverse flow velocity in this area has a decreasing trend and the values are quite close. In theory, this is a manifestation of flow diversion. Because of the existence of vegetation community, flow is blocked from advancing and flows towards the channel wall, resulting in a decrease in lateral flow velocity. Second, the simulation results and experiment results of downstream vegetation community do not fit well, and the simulation data need a long process to restore to the trend of the experiment data. Based on Fluent analysis, it can be seen that the lattice cannot well simulate the internal conditions of patch vegetation and the vegetation is too aggregated. During the simulation process, the flow cannot smoothly pass through the patch vegetation, resulting in a small flow velocity behind the patch and a large dead flow zone. However, the accuracy of Fluent simulation results in stable region is lower than that of LBM. Further improvement of lattice partition accuracy can ensure more accurate LBM simulation results, but will lead to multiple increase in calculation time.
In general, the simulated and experimental values are in agreement. This LBM is valid and reliable, and it can adequately simulate the flow of ecological open channels.

Parameter selection in ecological open channel LBM simulation
How to link the actuals to the simulations is another important issue in the LBM simulation. For the open channel flow problems, researchers believe that the simulated lattice Reynolds number and macro Reynolds number should be the same to reflect the actual flow problems (Mohamad 2011). Therefore, the parameters need to be adjusted in the simulation to achieve the same lattice Reynolds number as the macro Reynolds number. The macro Reynolds number can be expressed as where Re 1 is the macro Reynolds number, u is the inlet flow velocity of open channel, y is the width of open channel, and v is the kinematic viscosity of flow.
The lattice Reynolds number can be expressed as where Re 2 is the lattice Reynolds number, u LBM is the inlet flow velocity of LBM model, ly is the number of lattices in the y direction of the model, and v LBM is the kinematic viscosity of the model. The main parameters involved in the LBM are lattice size (∆), lattice fluid viscosity ( v LBM ), time step ( t ), relaxation frequency ( ), and fluid sound velocity (C). One of the fluid viscosities can be considered for the large eddy simulation method, and the specific calculation method can be found in Gac (Gac 2015). When using the LBM, only a choice of suitable parameters can ensure that the simulation matches reality and that the simulation results converge. When building LBM model, in order to keep LBM simulation consistent with macro problems, the following steps are adopted (Fig. 16).
• Calculate the macro Reynolds number according to Eq. (21). • Set the lattice inlet flow velocity u LBM (generally the order of magnitude is 0.1 or 0.2 to maintain the stability of the model). (22) • Preset the lattice size ∆ according to the macro problem to obtain the number of lattices in the y direction. • Ensure that the lattice Reynolds number Re 2 is equal to the macro Reynolds number Re 1 , and calculate the lattice fluid viscosity v LBM in combination with Eqs. (21) and (22). • Determine the fluid sound velocity C according to the macro problem and calculate the time step t ( t = Δ∕C). • Calculate the relaxation coefficient by Eq. (12) to determine whether the model converges and to adjust the mesh length (according to the simulation experience, the relaxation coefficient must reach 0.6 or more to keep the model stable).
According to the above steps, it can be found that the fundamental factor affecting the convergence of the model is the lattice size. Taking the dispersed vegetation open channel with a flow of 100 m 3 /h as an example, Fig. 17 shows the flow under different lattice sizes and relaxation coefficients. In Fig. 17a, the lattice length is 0.02 m, and the relaxation coefficient is 0.51. The lattice length in Fig. 17b is 0.005 m, and the relaxation coefficient is 0.8. The flow under different iteration times is shown in  Fig. 17a is four times that in Fig. 17b, so the number of iterations in Fig. 17b is four times that in Fig. 17a. When t = 50 in Fig. 17a, it is clearly seen that the flow is discontinuous and the change of flow velocity has a wavy shape. The flow is continuous when t = 200 in Fig. 17b. When the flow reaches the vegetation, it can be seen that the flow distribution in the simulation with a lattice length of 0.02 m is still discontinuous. The velocity has a snowflake-like distribution trend. When the grid length is 0.005 m, the flow is continuous and unobstructed. It can be clearly seen from Fig. 17 that lattice densification can make the simulated flow more realistic and the change of flow velocity can be more clearly displayed. When the relaxation coefficient is too small, the simulation results will even appear infinite or infinitesimal, and some values will appear in sporadic parts. Therefore, in the simulation, select the appropriate lattice size and determine the relaxation coefficient according to the macro problem, so as to achieve the required accuracy.
The lattice size should also be selected appropriately. If the lattice is large, the calculation speed will be accelerated, but the model may not converge. If the lattice is small, the model converges, but the calculation speed is very slow. Computing speed is often not proportional to the size of the lattice. Double the encryption of the model and the calculation speed is often twice or even more than twice as slow.
In general, the lattice Boltzmann equation is linear. In fact, the nonlinearity is embedded on the left side of the lattice Boltzmann equation. Therefore, LBM does not need to solve the Poisson equation at each time step to meet the continuity equation. This greatly reduces the calculation time. LBM can be regarded as an explicit method without simultaneous equations at each time step. Therefore, LBM has great development prospects in flow problems with ecological open channels.

Simulation of pollutant diffusion in compound channel
The research of Huang et al. (Huang et al. 2002) showed that the flood plain covered with vegetation will accelerate the flow velocity in the main channel. At this time, a turbulent mixing layer will be formed between the main channel area and the vegetation area due to the huge velocity gradient, thus greatly changing the river flow structure, and then affect the material transport characteristics in the river and change the mixed transport process of pollutants. For wide and shallow natural rivers, especially when the vegetation on the flood plain is submerged vegetation, the vertical and horizontal unevenness of water flow is relatively small. Therefore, the vertical unevenness of the flow can be neglected. In the vegetation area, the flow velocity is small due to vegetation resistance, while in the main channel area, the flow velocity is large. A strong mixing zone will be formed between the two, forming a strong mixing of material and energy in the flow body. Therefore, this section uses the experiment data in the compound open channel for verification. Firstly, LBM is used to simulate the compound open channel to obtain the basic flow information. Figure 18 shows  Figures 19 and 20 show the relative error between the measured value and the simulated value of the flow velocity in the compound channel with vegetation. The change trend of the simulated value is the same as that of the measured value. The average relative error between the simulated value and the measured value is 6.99%, and consistency is acceptable. Through the above analysis, it can be proved that the flow field information of the compound open channel with vegetation obtained by LBM simulation is accurate. A subsequent simulation verification can be conducted using it.
The results of Fig. 18 are input into the CA. The cellular size of the CA is the same as that of the LBM, so each cellular in the CA has its own flow velocity u t i,j and corresponds to each other one by one. Take the maximum flow velocity of the whole flow as u max . The flow influence coefficient of each cell can be calculated according to Eq. (5). Figure 21 shows the simulated and measured values at the upstream and downstream ( x = 2.5 and x = 5.5 ) of the open channel when y = 0 . The simulated value change trend is basically consistent with the measured value, and the fitting is reasonable. It indicates that the CA model can be used for the pollutant model of the vegetation open channel. The pollutants pass through the upstream measurement point in about 18 s and then move in about 25 s to reach the downstream measurement point. The pollution concentration of the upstream measurement point is significantly higher than that of the downstream, indicating that the pollutants are diffused in the horizontal direction. This is consistent with the meaning of diffusion coefficient in Eq. (4). According to this phenomenon, the movement of pollutants can also be described as the movement of still diffusion and superposition with water. Figure 22 shows the simulated and measured values of the upstream and downstream ( x = 2.5 and x = 5.5 ) of the open channel when y = 0.25 , where is the vegetation accumulation area. By comparing the simulated values with the measured values, it can be found that the CA can also better simulate the change of pollutant diffusion with time in the vegetation area. Figures 23 and 24 show the error between the simulated and measured values simulated using the CA. It can be seen in Figs. 23 and 24 that most errors between the simulated and measured values are within the allowable range. It shows that using CA to simulate pollutant transport is reliable. Some simulated values differ greatly from measurements, and the diffusion of pollutants is greatly affected by water flow, resulting in large random errors. In the experiment, vegetation is impacted by water flow and sways with the water, resulting in changes in the diffusion direction and speed of pollutants.
From Figs. 21 and 22, it can be found that the pollutants pass through the vegetation area longer. Most of the pollutants in the nonvegetation area pass through about 20 s, but the pollutants in the vegetation area pass through the area after about 40 s. It indicates that the presence of vegetation will lead to the accumulation of pollutants and slow down the diffusion velocity. Nepf et al. (2007) also found a similar situation. They considered that the flow blocking effect of vegetation affected the diffusion process of pollutants, because the presence of vegetation canopy would greatly slow down the mixing velocity of pollutants with flow. The vortex street behind the vegetation will capture pollutants, leading to the accumulation of pollutants. The amount of pollutants captured by the vortex is related to the Reynolds number of the flow (White and Nepf 2003). Gerrard (1978) pointed out that in the absence of vortex shedding, the capture and release of this region occurred through diffusion. The average residence time of pollutants in this region is t ∼ d 2 ∕(4K) , where K is the molecular diffusion constant. When there is vegetation in the open channel, eddy current will be generated, and the eddy current will fall off with a frequency of f s . This phenomenon controls the capture and release of pollutants in the vegetation wake area. So the average residence time of pollutants in this area is t ∼ 1∕f s . In both cases, the residence time of pollutants is longer than the advection time, t > d∕u . This conclusion is consistent with the simulation results. The time of pollutants passing through vegetation areas is longer than that passing through nonvegetation areas.
Compared with other pollutant diffusion models that need to solve the partial differential equations simultaneously, this model can realize the efficient simulation of pollutant diffusion in vegetation channels. The LBM-CA model has high computational efficiency and simple evolution rules. The standard limits of specific pollutants can be embedded in the LBM-CA model. During the simulation, the number of cellular with heavy pollution, light pollution, or trace pollution in the river at a specific time point can be read, providing solutions for the accurate treatment of river pollution.

Conclusion
Cellular automata and lattice Boltzmann method are combined in this paper to construct a pollutant diffusion model  Data availability Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Declarations
Ethical approval Not applicable.
Consent to participate Written informed consent for participation was obtained from all participants.
Consent for publication Written informed consent for publication was obtained from all participants.

Competing interests
The authors declare no competing interests.