Backward Integration of Diffusion Equation

When the parabolic differential equation is integrated backward in time, it can create unwanted shortwaves with large amplitude. Hence, instead of solving it as a differential equation, the diffusion equation is converted to the equations of volume-integrated-concentration, and mixing/diffusion is treated as subgrid-turbulent-fluxes across the cell boundaries. Those equations become a set of linear algebra equations and can be solved in both forwardand backward-in time. The proposed method has been validated by the numerical simulations of an idealized case, which consists of 5 different sizes of concentric cylinders with different species. The time evolution of compositions shows that the concentrations in each cylinder can change drastically with time. For the data collected at downwind region, the proposed reverse-in-time integration can be used to assess the concentrations at the source regions, which can be quite different from those derived from the conventional backward-trajectory method without mixing. It also shows that the traditional forward-trajectory or backward-trajectory method without mixing (i.e., Lagrangian method) widely used in meteorology and air pollution can misinterpret the property of fluid parcel at both upwind and downwind regions significantly. Keyword: Forward and backward integration; Diffusion; Parabolic and hyperbolic equations; Mixing/Diffusion-time scale; Turbulence; Pollution.


INTRODUCTION
Turbulent mixing and diffusion are important in fluid dynamics.It has been included in most atmospheric and ocean models to simulate/predict weather, climate, and oceanic phenomena.Numerical models are also applied to study the transport and dispersion of aerosols and pollutants in the atmosphere and oceans.The advection-diffusion equations in those models are usually integrated forwardin-time, starting from the source as an initial condition to simulate the transport and dispersion of the plumes or thermal as time increases (Sun, 1989;Sun et al., 2013a, b).From the concentrations in the downwind region, backwardtrajectory-method (Choi et al., 2001) is frequently used to estimate the property of the parcel at source derived from backward advection equations without mixing, because the reverse diffusion can easily become unstable for short-wave perturbations.It is also noted that Cauchy problem for the reverse heat (diffusion) equation (∂u/∂t = a∂ 2 u/∂ 2 x, wih a < 0) is incorrect in the sense that it does not have continuous dependence of the solution on the initial function (Lavrentiev et al., 1980).Small errors of the data lead to a large errors in the results.It is necessary to filter out all information containing higher Fourier coefficients (Engl et al., 2000).Hence, Scientists have proposed "Regularization" by introducing addition information in order to solve this ill-posed problem.For example, Sikora and Palka (1980) proposed a method consisting of solving the Volterra integral equation of the first kind by means of Tikhonov regularization method.This approach requires complicated mathematics and computation, which may not be applicable to the real pollution problems.
It is also noted that the standard linear parabolic models based on Fourier's law (1822) or Fick's law (1855) predict an infinite speed of propagation, which may violate physics.Hence, Cattaneo (1958), Vernotte (1958), Chester (1963), and others proposed that Fourier equation should be upgraded from the parabolic to a hyperbolic form, where α is called the speed of second sound.Then, the equation can be solved as waves, which can propagate forward and backward.However, the solutions are sensitive to the relaxation time (τ r = κ/α 2 ) and may become complex numbers.The solution may be acceptable according to the second law of thermodynamics, it seems unrealistic for gaseous material and possibly for most solids (Taiel, 1971).Hence, forward-in-time integration was applied in most traditional diffusion equations (Sun, 1982).
As shown in Sun and Chang (1986) and Sun (1988Sun ( , 1989)), the change of concentration C can be rewritten as: where V is the velocity, and C'V ' is the subgrid-turbulent flux of C. We can integrate (1b) over any specific cell if the turbulent flux around the boundaries can be formulated properly.Hence, we can convert (1a) from a differential equation into equations of volume-integrated-concentration, which can be integrated forward or backward in time.The proposed method can avoid creating undesirable shortwaves, because the size of the cells are well defined.If there is no internal source or sink, according to conservation of mass, the change of mass in each cell is equal to the net mass fluxes across the boundaries, consisting of cell-and subgrid-scale fluxes.The diffusion process can be interpreted as the subgridflux, which depends on density perturbations and the turbulent kinetic energy (TKE) at the boundaries.For simplicity, it is assumed that size of each cell is time independent and density inside each cell is uniformly distributed, then, we can easily convert the integration equations to a set of linear algebra equations, which depend on the size of cells, concentration, and the turbulent velocity.The method can be modified and applied to the cases that containers, wind, and turbulences are function of time.The numerical results for the idealized cases are quite accurate if integration is comparable or less than the diffusion-time-scale, τ, which is a function of volume (V) and TKE.If integration is longer than diffusion-time, majority of the original material in the individual cell is replaced by the substances of the surroundings in the forwardin-time integration.It is also noted that the error increases as the reverse-integration time increases.
Severe air pollution usually occurs when the atmosphere is very stable with little mixing, under such condition we can obtain highly accurate results after a long backward-in-time integration.If the atmospheric stratification is near neutral with strong mixing, pollutants can disperse quickly, and air quality may become more tolerable, under this condition we cannot perform a long reverse-integration without introducing significant errors, because the time scale of backwardintegration is limited by the mixing/diffusion-time-scale.Since the property of the air parcel at the downwind regions can be significantly different from that at the source regions.Hence, the traditional forward-backward trajectory method without including mixing is valid only when the time scale is much shorter than the diffusion-time-scale τ, which will be defined later.

EQUATION AND NUMERICAL MODEL
According to Leibniz' rule, the change of density j m  = ' (cell-sized mean)+ (perturbation) for the m th specie in the j th cell with velocity V j = V (cell-sized mean)+ j V '(perturbation) where is the sink of j m  , it is also assumed that V j  = 0 and sizes of the cells remain constant.The bar of cell-sized mean will be dropped hereafter.Following the Lagrangian frame, the change of the mass  for the m th specie in the j th cell with time- independent volume V j can be given by: , , , , ' ' 1 for i = 1, 2, .., J; j = 1, 2, .., J; and m = 1, 2, .., M where A i,j is the surface area between the i th and j th cells, and , ' i j u is the turbulent velocity on A i,j .The terms on the right hand side indicate the material of m th specie getting out from j th cell, getting into j th cell from outside, and the internal sink inside j th cell, respectively.For simplicity, V j , A i,j , and , ' i j u are time-independent in this study.Otherwise, Eq. ( 2) should be applied.The finite difference form of Eq. ( 3) becomes for j = 1, 2, .., J; i = 1, 2, .., J; and m = 1, 2, .., M An implicit-time step is applied on the right hand, where q = 1 -p, and 0 ≤ p ≤ 1. (4) Becomes: We assume that p = 0.5 in order to have the 2 nd order accuracy in time.Eq. ( 5) can be solved by LU-method for the forward-in-time from an initial condition to study the dispersion of species as time increases.It can also be solved for the reverse-in-time from the concentrations collected from the downwind regions to assess the initial conditions and time evolution between the time intervals.
The model consists of two layers of cylinders as shown in Fig. 1.The mean height of the PBL is somewhere between 1 km to 2 km (Medeiros et al., 2005), but it may be higher over the industrial and urban areas, so we set it as 2000 m for the depth of the lower layer, which is topped by the cylinder E (light yellow) with radius of 30,000 m and thickness of 2,222.2m, in the lower part of the free atmosphere above the PBL.The lower layer consists of 4 concentric cylinders: A (black), B (red), C (green), and D (blue), with radii of R AB (between cylinders A and B) = 5,000 m, R BC (between cylinders B and C) = 10,000 m, R CD (between cylinders C and D) = 20,000 m, and R D (outside of D) = 30,000 m.The volume is 157,079 × 10 6 m 3 for cylinder A; 471,239 × 10 6 m 3 for cylinder B; 188,496 × 10 7 m 3 for cylinder C; 314,159 × 10 7 m 3 for cylinder D; and 628,256 × 10 7 m 3 for cylinder E. Different species with uniform density are assigned to different cylinders initially.The turbulent velocity u' between two neighbor volumes should be proportional to the TKE, which can be obtained from observations or numerical models but is not calculated in this study.Hence, an assigned value of TKE is used here.It is noted that the TKE was explicitly calculated in many atmospheric models, including Purdue Regional Climate Model (Sun and Chern, 1993;Sun et al., 2009;Min and Sun, 2015) and NTU-Purdue Model (Hsu andSun, 2001, Sun andSun, 2015), and CRESS model (Tsuboki and Sakakibara, 2007;Sun et al., 2012Sun et al., , 2013c)), etc.
Severe pollutions usually occur in a calm, stable condition.Observations showed that for the weak-wind (< 5 m s -1 in the lowest 200 m), very stable boundary layer (with Richardson Number ≥ 0.6).The boundary layer was very shallow (sometime < 10 m deep), and turbulent fluxes between the earth's surface and the atmosphere were found to be essentially shut down with TKE being near zero (Banta, 2008).The PBL above this strong surface inversion may be still quite stable, hence, the value of horizontal u' is set to 2.5 × 10 -2 m s -1 .The vertical mixing u' between the cylinder E in the upper layer and those in the lower layer is set equal to 0.1 of the horizontal u' in the lower layer, because vertical mixing is suppressed under the stable stratification.Since step functions are applied to j m  , density perturbation across the cells' boundary, '

Forward-in-Time Integration
In order to simplify tracing the dispersion of the material with time, each cylinder has its unique specie with one unit amount of density initially, as shown in Table 1.It is also assumed that the specie will maintain its identity during integration.The time evolutions of the density of species in each cylinder are shown in Figs.2(a The concentration of 1 A   0.2 at t = 2 × 10 5 s shows that most of original material is replaced by the substances from outside by this time, especially those with vast volumes, only 7% of the original material remains at t = 4 × 10 5 s.
On the other hand, the concentration of A at tf in the downwind will be interpreted as that of A initially according to the backward-trajectory method.Because of diffusion, the material of each cell in the downwind is quite different from the original source.Hence in meteorology and other fields, , (where H L is the depth of lower cylinders), is same for each lower cylinder, the increase of 5 5 5 5 , , , and are also the same, although the mass increment of each cylinder depends on the volume of cylinder.

Reverse-in-Time Integration
Using the density distribution of each specie at t = 4 × 10 5 s shown in Table 2 as initial condition, we integrate Eq. ( 5) reverse-in-time, i.e., we calculate ( )   by Eq. ( 5).The integrated density from t = tf = 4 × 10 5 s to the end at t = ti = 0 are shown in Figs.3(a)-3(e) for species 1 to 5, which are indicated by black, red, green, blue, and light blue lines marked with "○" signs.
Figs. 3(a)-3(e) show that the reverse-in-time integrations (with ○) coincide the track of the forward-in-time integrations (with +).The results end up almost identical to the initial condition used in the forward-in-time integration.ρ 5 A (Fig. 3(e)) has an error of 1.15% after 8,000 times of forward-Fig.3. Forward-integrations + and backward-integrations ○ for species among different cells: (a) 1 1 1 1 1 , , , , and ; , , , , and ;   integration can handle the counter-gradient flux, which is difficult to simulate in the diffusion equation even in the forward-in-time calculation (Sun and Chang, 1986;Sun, 1988Sun, , 1989)).
The simulations show the importance of mixing among the cells.Because of mixing, the material at the downwind region can be significantly different from its original source, which cannot be derived from the forward-trajectory method without mixing.Meanwhile, it is impossible to predict the property at the sources from the downwind regions using the conventional backward-trajectory-method.Both forwardtrajectory and backward-trajectory methods without mixing are referred as the Lagrangian method and have been widely used in meteorology and air quality researches.
Here, we provide a different approach.From the observed concentrations at downwind and meteorological conditions, the proposed backward-in-time integration method can be applied to analyze the contributions coming from the source and surroundings in the upstream regions.

Time Scale of Mixing/Diffusion
Both forward-and reverse-in-time integrations, the solutions depend on where L j is the length scale of the j th -cell, τ j , (   ) is the time scale of mixing/diffusion for j-cell, and n is the number of time-step of integration.The time scale of mixing/diffusion τ j ~ L j /u' ~ 200,000 s for Cylinder A, and τ j ~ 600,000 s for cylinder D. Table 3 shows that the error combining both forward-and backwardintegrations decreases with the increase of cylinder size, as expected.The error for cylinder A is about 0.01, and 0.001 for cylinder C, both errors are associated with specie 5.The results are quite accurate if the integration time is order of τ j or smaller.However, if the integration time is longer, the error may increase drastically.
If the size of cylinders remains constant, we can obtain the same results as long as , ' i j n u t does not change.Hence, the integration time t n should be shorter if the turbulent velocity , ' i j u is stronger in order to keep t n /τ j ≤ 1.On the other hand, t n can be longer for weaker turbulence.For example, if the volume of cells remains the same.The results shown in Table 1, 2 and Figs. 2, 3 for the case with tf = 4 × 10 5 s and 1,2 ' u = 2.5 × 10 -2 m s -1 are also applied to the cases with tf = 2 × 10 5 s and 1,2 ' u = 5 × 10 -2 m s -1 ; or tf = 6 × 10 5 s and 1,2 ' u = 1.67 × 10 -2 m s -1 , etc.A larger value of , ' i j u implies stronger mixing among the cylinders.The property of cylinders dilutes quicker, especially for the small-sized cylinders.On the other hand, a smaller value of , ' i j u corresponds to weaker mixing, such that original elements can retain longer.
If the turbulent velocity keeps constant, the results remain the same when    remains constant according to Eq. ( 6).It indicates with the same turbulent intensity, a small sized cell quickly mixes with surrounding and is more difficult to trace back from reverse integration.Hence, the results shown in Figs.3(a)-3(e) and Table 3 are also valid for the cylinders with different sizes or different depth of PBL, as long as time-integration is adjusted to length scale accordingly (i.e.,   . Hence, a longtime integration should be applied to the large-sized cells under weak turbulent condition.Since each specie remains its identity (i.e., without chemical reaction), the results obtained here can be applied to the case such that each cell has any different combinations of those species initially.
It is noted that wind, turbulent velocity and the size of cells in the real atmosphere and oceanography are timedependent.Hence, the turbulent velocity and the volume of cells proposed here should be modified according to the atmospheric model simulations or observations.Similarly, sink/source term should include deposit and chemical reactions in order to simulate the transport and dispersion in the real atmosphere and oceans, as discussed in Sun et al. (2013a, b), andYang (2004).

CONCLUSIONS
Backward integration of the parabolic equation can easily create undesirable large amplitude of the short waves.
Hence, equations of volume-integration of concentration are proposed, in which diffusion is treated as the subgridturbulent-fluxes across the boundaries of the cells, which depends on the turbulent velocity and density perturbation.For given turbulent velocities, the change of density of different species in each cell can be converted to a set of linear algebra equations and solved for both forward-and reverse-in time.The proposed method has been validated by the numerical simulations, which consist of 5 different sizes of concentric cylinders with different species initially.The results reveal the time evolution of compositions in each cylinder.Simulations also show that the material in each cell can change significantly during integration.The reverse integration goes through the counter-gradient diffusion process to reproduce the initial concentration of each cell, which can be very different from the downwind regions.Since the conventional forward-or backward-trajectory method ignores mixing process, it cannot be used to predict the property at the source from the downwind region, and vice versa.The method proposed here can also be applied to different cell-size, or different turbulent velocity, as long as 1/3 , ' / ( ) i j n j u t V is adjusted to the real environments.The results also show that the error increases as the reverse integration time increases.However, if integration time is comparable or shorter than the mixing/diffusion-time scale, the method is quite accurate.
than that with a continuous-variation of density among different cells.Since flux at the boundary is determined by ' ' u  , a smaller value of u' may compensate a larger value of ρ' applied in the study.Different value of u' will also be discussed later.The time interval Δt = 50 s is applied to both the forward and reverse integrations.It is noted this method proposed here can be modified and applied to the cases when u' and the size of cylinders change with time according to (2).

Fig. 1 .
Fig. 1.The size and arrangement of cylinders in this study.


)-2(e).The black color indicates the material originally in cylinder A; red in cylinder B; green in cylinder C; blue in cylinder D; and light blue in the upper cylinder E. Initially, 1 A  (specie 1 in cylinder A, black line) decreases with time exponentially, which is replaced by the material coming from cylinders B in light blue).The density of those invaders starts from zero and increases with time as shown in Fig. 2(a).The decrease of 1 A  slows down as time increases, because amount of ρ 1A decreases with time.32 as shown in Table1and Fig.2(a).
t = 0.It indicates that the reverse-

Table 1 .
Initial value of k m .

Table 3 .
Value of ρ