Analytical and Numerical Solutions of Radially Symmetric Aquifer Thermal Energy Storage Problems

The paper discusses analytical and numerical radial solutions of the differential equations for heat transport in water-saturated porous media. In particular, a similarity solution is obtained for a 2D-horizontal confined aquifer with constant radial flow. Numerical solutions are derived using a high-resolution Lagrangian approach suppressing spurious oscillations and artificial dispersion. There is a good correspondence between numerical and analytical solutions. The primary purpose of the investigation has been to calculate the recovery factor of an Aquifer Thermal Energy Storage 5 (ATES) system with a cyclic repetition of injection and pumping. Solutions covering both instantaneous and delayed heat transfer between fluid and solid, as well as time varying water flow, are derived and applied to a one-well test case. In hydrological terms, these solutions are relevant for a wide range of problems where groundwater reservoirs are utilized for extraction and storage (viz. irrigation; water supply; geothermal extraction).


Introduction
Heat transfer in porous media has received considerable attention and is the topic of a number of investigations during the last decade (Bejan and Kraus, 2003).A driving force for research on this subject is engineering applications, such as geothermal systems (Ganguly and Kumar, 2014), heat exchangers (Diao et al., 2004), thermal insulation (Birhanu et al., 2015a;Kim et al., 2010), and safety issues regarding storage of nuclear waste (Sun et al., 2010).
In addition to equations for the fluid flow, the mathematical model of heat transfer in porous media is given by second-order partial differential equations for heat energy conservation and flow in the model domain.Two kinds of models can be applied to investigate thermal characteristics of conduction and advection within a porous medium, namely, a thermal equilibrium model and a thermal nonequilibrium model (Yang et al., 2011).The difference between the two models is the thermal coupling between the liquid and the solid phase.For the equilibrium model the coupling is modeled as an instantaneous heat transfer.This assumption is close to the reality for homogeneous aquifers with solid particles of minor size (diameter d p < 1 mm).For the non-equilibrium model there is a time delay attached to the heat transfer between the two phases.In the literature of transport in porous media, this model is usually called a double porosity model which may also be expanded to a dual permeability classical dipole geometry.In this case, the flow field may be simplified to axis symmetry where the flow velocities are governed by the injection/pumping rate and the aquifer porosity.In the present study, we take advantage of a simplified flow field and solve the transport equation by simple analytical and numerical methods to evaluate the energy efficiency of an idealized ATES system.Numerical solutions of transport problems are usually affected by artifacts, and because all mathematical models are simplifications of reality, boundary conditions should be specified with great precaution.Birhanu et al. (2015b) showed that the most evident boundary condition of temperature at the top of an unconfined aquifer gave unphysical energy efficiency for simulation experiments of a real ATES system.In this study we use analytical solutions to understand the quality of numerical simulations by doing (simple) numerical experiments.One such experiment is an idealized ATES production sequence consisting of repeated injection and pumping of hot water in a confined aquifer.The performance of the alternative solutions is quantified by a recovery (or efficiency) factor.The results are presented in Section 5.

Mechanisms and Equations
A thermo-hydraulic analysis requires calculation of simultaneous water and heat transport in an aquifer consisting of a solid porous medium (s) with pores filled with water (w).The water flow depends on properties of the water as well as the solid, and the gradient of the hydraulic head, as stated in Darcy's law (Molson et al., 1992;Nield and Bejan, 2013): Here, q is the specific discharge or the Darcy velocity, k the intrinsic permeability tensor, z the elevation of the piezometric head relative to a datum level, p the fluid's pressure, ρ water mass density, g the acceleration of gravity, and µ the dynamic viscosity of water.Furthermore, K = ρgk/µ is the hydraulic conductivity, and φ = p/ρg + z is the hydraulic head.The volume average velocity differs from the velocity of the water in the pores, the so-called seepage velocity, u = q/n, where n is the (effective) porosity (Kangas and Lund, 1994).The water density and, more pronounced, viscosity vary with temperature.However, we shall here assume that the flow field represented by q remains independent of the temperature changes.Assuming that the solid matrix as well as the water are incompressible, mass conservation combined with Darcy's law leads to the Poisson equation for the hydraulic head, where Q w (t) is a source/sink term.
The heat energy content per aquifer volume unit may be written where c is the specific heat, and subscripts w and s refer to water and solid.At a local temperature equilibrium where, T w = T s = T , the heat content may be expressed as (ρc) m T , where see Kangas and Lund (1994) and Hecht-Méndez et al. (2010).In the following we will use the convention ρ m c m for (ρc) m .
The water flow causes advection of heat, whereas conduction/diffusion of heat takes place both in the solid and the liquid, and λ w,s are the heat diffusion coefficients.If the two media are at a local thermal equilibrium, the volume average diffusive heat flux may be expressed by where λ m is a bulk aquifer heat diffusion coefficient, see Kangas and Lund (1994) and Nield and Bejan (2013).Other expressions for λ m , e.g.porosity-weighted geometric and harmonic means are also discussed in the literature (Nield and Bejan, 2013).In addition, the heterogeneity of the pores induces a certain amount of thermal dispersion, parametrized, in its simplest form as Here α is the thermal dispersivity length, and the total diffusion flux becomes q T + q d (Bakr et al., 2013).
If a thermally insulated aquifer with initial temperature T w = T s and q = 0, energy conservation implies that the system attains an equilibrium temperature T m equal to the weighted mean How fast the temperature equilibrium between water and solids is reached depends on the efficiency of the energy exchange between the two media.It turns out to be reasonable to express the heat exchange per time and volume unit as where h is a heat transfer coefficient (Nield and Bejan, 2013;Kreith et al., 2010).The coefficient varies with temperature and flow, in particular for large flows.Following the discussion in Nield and Bejan (2013), h may be expressed as h = a ws h v , where a ws = 6 (1 − n) /d p is the surface area of the water/solid interface per volume unit.For low Reynolds numbers, h v may be expressed as Here, d p is the size of the grains making up the solid.The expression for h then becomes A rough estimate of the time scale ∆t towards thermal equilibrium may be obtained from the energy exchange per time unit at the start of the heating, P = h (T w − T s ), compared to the required amount of energy to be transferred, The time scale is thus only dependent on basic material constants and the grain size.With typical values for rock, we obtain A similar time scale may actually be derived from the heating of spheres discussed in Gockenbach and Schmidtke (2009).
For an elemental aquifer volume R with boundary ∂R, the solid's integral conservation law reads Note that since the solid is not moving in this case, there is no advective heat for the solid phase.Similarly, the integral conservation form for heat in the water is The differential forms of the conservation laws with the assumptions above become We observe that when T w is kept constant and diffusion is neglected, the natural time scale (inverse rate of change) in Eq. ( 19) is essentially ∆t.For d p less than about a millimetre, the thermal equilibrium is virtually spontaneous and we may assume that T w and T s are equal.
For the case where T = T w = T s , we obtain by adding Eq. (17-18), and the corresponding differential form, If the parameters like c, ρ, λ and the flow q are assumed to be independent of T , then dividing through with c m ρ m in Eq. ( 22) leads to In our model, the flow q is caused by water injected or pumped with a discharge rate Q(t) from a well located at the origin.
When Q(t) > 0, water is injected from the well into the aquifer, causing a flow away from the well.During pumping, Q(t) < 0 and the flow is directed towards the well.Utilizing symmetric geometry of the aquifer near the well, the flow is q = q d i r in which the discharge velocity The width W and height H are constants characteristic for the aquifer.In this case, Eq. ( 23) becomes Using the same symmetry considerations in the nonequilibrium case Eq.(19-20) become temperature) on the other hand, are usually more easy to monitor.In such cases, the physical parameters are estimated by solving the inverse problem.Before the advent of computer technology, this was done by dimensionless solutions of the analytical expression, which provided tables or so-called type-curves.The inverse problem was solved by curve fitting of the empirical data to the analytical type curve.Today, analytical solutions are applied in similar ways, but the curve fitting is substituted by numerical perturbation of the involved variables.Sensitivity analysis is another example where application of analytical solutions is convenient.After estimation of optimal parameters, the relative impact of the uncertainties might be evaluated by simple numerical perturbation of the involved parameters.Here, in this context, the motivation for using analytical expressions of the transport equation is due to the numerical challenge of solving the transport equation if advective flow is dominant.In such cases, numerical solutions are prone to numerical dispersion.Even though analytical solutions simplify real transport, numerical artifacts do not affect the solutions.Therefore, by using the same simplified velocity field for both the analytical and the numerical solutions, the performance of the numerical algorithm can be evaluated directly by using the analytical solutions as benchmarks.
We shall consider the formation of a hot water plume in a local thermal equilibrium aquifer generated by a constant hot water source at the origin.Consider Eq. ( 25) for convenience normalized such that κ d = 1 and with initial and boundary conditions, If the diffusion term is negligible, Eq. ( 25) becomes a simple hyperbolic equation which, for any initial temperature distribution In particular, for the boundary conditions in Eq. ( 28), the hyperbolic solution is the moving front T (r, t) = H c r d − dκ d t where H c (x) is the complementary Heaviside function (= 1 for x ≤ 0, = 0 for x > 0).Since ρ w c w is typically about twice as large as ρ s c s , the ratio ρwcw ρ m c m depends on the porosity n and varies between 1 and 2. The temperature front thus moves significantly faster than the discharge velocity q, but more slowly than the average seepage velocity, u = q/n.
The 1-dimensional case When d = 1 Eq.( 25) becomes where we choose the more standard spatial variable x rather than r.In this case, a well-known similarity variable is η = (x − κ 1 t)/ √ λt, which, inserted into the equation results in with the general solution T (η) = C 1 erf(η/2) + C 2 .However, no solution from this collection satisfies the boundary condition T (0, t) = 1.Nevertheless, the similarity solution of the closely related problem satisfying the initial values T (x, 0) = H c (x), x ∈ R is a good approximation: A modification of this solution, satisfying all conditions in Eq. ( 28) exactly has been derived by Ogata and Banks (1961); see also Bear and Cheng (2008), Eq. 6.4.30: The similarity solution Eq. ( 31) and the exact solution Eq. ( 32) are presented in Fig. 1 together with the hyperbolic front solution T (x, t) = H c (x − κ 1 t).As expected, the similarity solution in Eq. ( 31) does not satisfy the boundary conditions at x = 0. Still, as λ tends to 0, the solution approaches the hyperbolic front, and Eq. ( 31) becomes a very good approximation.
The 2-dimensional radial symmetric case For a 2-dimensional problem, assuming radial symmetry, Eq. ( 25) becomes which may be rewritten as Again, it turns out that assuming the similarity variable η = r/ √ λt, we obtain an equation with general solution The solutions may be written in terms of the incomplete Γ-function, defined as The radial 2D similarity solution becomes or and the solution is shown is Fig. 2.
The spherical symmetry 3D-case is easily seen to have intrinsic scales for r and t involving κ 3 and λ, and no simple similarity solution exists.The scaled equation may however be transformed to a 1D heat equation with a space dependent diffusion coefficient (Philip, 1994).Actually, the numerical algorithm below applies a similar transformation.

Computational procedure
We shall now consider a numerical algorithm for solving Eq. (19-20), or the equilibrium model Eq. ( 23) in the case of symmetric geometry.It turns out that the numerical solution of these problems is nontrivial.They are typically advection dominated, and we have already seen in the previous section that the temperature profile is a sharp front moving away from the source.
In the radial and spherical case, the flow becomes very large close to the origin, leading to an almost hyperbolic equation in this region.Advection dominated problems are notoriously difficult to solve numerically.Popular schemes, like central differencing schemes result in unstable or spurious oscillatory solutions.Upwind discretization for the advection term avoids oscillations, but does create artificial diffusion, leading to a smoothed temperature front when applied on a coarse grid.Several other methods have been proposed and discussed in the literature, see e.g.Strikwerda (2004) and LeVeque (1992).We propose a computational procedure utilizing the special structure of the Eq. ( 25) or Eq.(26-27).The procedure is composed by the well-known Lagrangian approach combined with a coordinate transformation.The idea will first be explained for the equilibrium case Eq. ( 25).
Before discussing the numerical scheme, it is convenient to say something about the scaling of the problem.Consider a time scale T for the time t, and R for the space variable r.A reasonable relation between the two scales is R d = κd T , R is the distance a hyperbolic temperature front moves in time T for a constant κ d = κd .The temperature T will typically be scaled with the temperature of the injected water.With these scales we get the dimensionless equation where As a curiosity, notice that in the case of a constant κ d , by choosing R and T such that κ d T /R d = λT /R 2 = 1, the parameters are completely absorbed, in the sense that also λ becomes equal to one.This is possible for d = 1 and d = 3, but not for d = 2.However, we will not pursue this curious issue any further here.
To handle the singularities in origin in the radial and spherical cases, we introduce the transformation s = r d /d so that ∂s = r d−1 ∂r, valid for all dimensions d.In this case, Eq. ( 40) becomes Notice that a(0) = 0 for the d ≥ 2, elucidating the hyperbolic nature of the problem close to the origin.The numerical difficulties of hyperbolic problems can be resolved by using a Lagrangian method: Given a path s(t) in the (s, t) plane.The solution along this path is T (s(t), t) and the total derivative of T with respect to time becomes If we let the path s(t) satisfy ds/dt = v(t), the advection term is eliminated.In fact, the paths s(t) are the characteristics for the hyperbolic equation we obtain for λ = 0.As a result, Eq. ( 41) can be solved as a system of differential equations: The first equation is an ordinary differential equation, whereas the second one is a heat equation with a space dependent diffusion coefficient.This can be discretized in space by some appropriate finite difference schemes, e.g.
with a i+1/2 = (a(s i+1 ) + a(s i ))/2, and initial values s i (0) = s i,0 and T i (0) = T (s i,0 , 0).The procedure is significantly simplified if v is constant, in which case the characteristics s(t) are just straight lines.
The spacial domain can be extended to R by defining a(s) = 0 for s < 0. In this case, we may solve Eq. ( 41) with the boundary conditions: When water is injected, v > 0 and the temperature of the water at the well is T (0, t) = 1.This is realized by choosing T (s, t 0 ) = T 0 (s) whenever s > 0 and T (s, t 0 ) = 1 for s ≤ 0.Here t 0 is either the initial time or a switching time, that is whenever v(t) changes from negative to positive (from pumping to injection).The procedure is illustrated for the injection phase in Fig. 3.
In order to be able to resolve a sharp front, the characteristics s i (t) used in the discretization can be concentrated around it.
The initial computational domain is (−S int , S int ) where S int is chosen sufficiently large to avoid any influence from the boundaries.The concentration of characteristics around the front is achieved by using where p is a positive integer, the higher p the stronger concentration.In our experiments, we have used p = 3, N = 100 and S int = 2.4.The underlying system of ordinary differential equations is solved by a standard solver in MATLAB (ODE15s).
For comparisons, Eq. ( 40) is also solved by a standard difference scheme with constant stepsize.In this case the advection term is approximated with an upwind scheme, (∂T /∂r)(r i , t) ≈ (T i (t) − T i−1 (t))/∆r.For the diffusion term a central difference scheme is applied.The spatial gridsize is ∆r = 0.012.The problem is solved for different values of λ, and the results are shown in Fig. 4 together with the exact solutions given by Eq. ( 39).For λ = 0.1, the diffusion is large and the artificial diffusion of the upwind method is insignificant.For smaller values of λ, the front remains sharp, and the effect of the artificial diffusion of the upwind method becomes quite pronounced.
The Lagrangian approach preserves the sharp temperature front.

The non-equilibrium case
We now consider the non-equilibrium case Eq. (26-27), which in scaled form becomes where Again, for a given time scale T it is appropriate to choose a spatial scale R such that v(t) is of the size of 1.The boundary conditions in the injection case (v(t) > 0) are (with T inj = 1 if the temperature is scaled).
By applying the transformation s = r d /d and using the Lagrangian approach, we can find the solutions T s (s(t), t) and T w (s(t), t) on the characteristics s(t) from where u(t) is the velocity of the temperature front, typically u(t) = ρ w c w n/(ρc) m v(t).The formulation is used to construct a spatial grid which is dense around the steep solution profile, and moves with it.This is illustrated in the following example.
Example 2. Consider the nondimensional equations Eq. (48-49), using d = 2 and parameters and a time dependent flow, v(t) = cos(πt).For 0 < t < 1, both injection, 0 < t < 0.5, and pumping, 0.5 < t < 1, are demonstrated.The equations are solved with the Lagrangian approach, using a central difference approximation for the diffusion terms, and a downwind scheme for the advection term.We have used a spatial stepsize ∆s i , varying from 1.1 • 10 −6 to 0.14 and concentrated around the temperature front of the water.The semidiscretized system is solved by MATLABs ODE15s.
The result of the simulation is presented in Fig. 5.We can clearly see how the hot water plume develops with time t.The concentration of characteristics moves with the temperature front, which remains sharp.At the same time, there is a heat exchange from water to solid, and the temperature of the water behind the front is reduced at the same time as that of the solid temperature increases.The diffusion is almost negligible most of the time, but will cause a slight smoothing of the front.When t → 1, the front approaches the well with increasing velocity, and the smoothing of the front becomes more pronounced.This example illustrates the numerical challenge, namely the coupling of advective and diffusive physics.Usually, the numerical solutions of such problems are suffering from artifacts, such as numerical diffusion and/or oscillations, but in this case it is possible to suppress the numerical artifacts to negligible levels.
5 Case study: Temperature profile near a well in an ATES system In this section, we consider the temperature propagation around a hot well in an ATES system.The physical and thermal properties of the Gardermoen aquifer were obtained from Goshu and Omre (2003) and Tuttle (1997), and is given in Table 1.The ATES system is typically operating in one of two modi: Injecting hot water during daytime and extracting it during nighttime, resulting in a full operational cycle of 24 hours.Alternatively, the warm water is injected during the summer months, and extracted in the winter, giving a cycle of one year.In reality, a combination of these is used, but in our study, we only consider the first modus, assuming the injection and extraction periods to be of equal lengths, 12 hours.
The Gardermoen aquifer is a delta structure deposited in a glacio-fluvial/glacio-marine environment during the last deglaciation of the Scandinavian crust (approx.10.000 B.P., Tuttle (1997)).The river discharge and the sediment load from the melting glacier were significant, which explains the wide range of grain sizes of the sediments from boulders (d p > 500mm) at the proximal side of the delta, to fine sand and silt (d p 1mm) at the distal side of the delta.The ATES system for this case study was located in the delta foresets with homogeneous fine sand, but it is interesting to compare the energy efficiency of this ATES with alternative locations.We therefore let the sediments vary form d p = 1mm, which corresponds to a foreset location, to d p = 500mm mimicking a location close to the glacial portals.In that case, the aquifer permeability would have been better, but to keep the experiment as simple as possible, we let the pumping rate and the porosity be the same for all grain sizes.In Table 2 the values of the heat transfer coefficient h, Eq. ( 14), and the time scale towards thermal equilibrium ∆t, Eq. ( 16), for different particle size are shown.So we can conclude that within time scales given by the injection/extraction periods, there is almost thermal equilibrium for realistic particle sizes.It is still of interest to see what happens in the initial injection phase, before thermal equilibrium is established, so we will solve Eq. (26-27).In a horizontal confined aquifer we can assume radial symmetry in the vicinity of a well, so d = 2. Initial and boundary conditions for the first injection phase is The equations are solved by the numerical approach outlined in Section 4.

Transient injection phase
In Fig. 6 we present the temperature profiles for the first few seconds of the injection period.The first row shows the situation for particle size d p = 1mm.Thermal equilibrium happens almost immediately in this case, but the energy transfer still has an effect in the sense that the temperature front become smoother.Also notice that after 0.15 sec, the solid temperature at the wall has reached to about 2/3 of the water temperature, while the water is almost cooled at the front.This is consistent with the fact that (1 − n)ρ s c s /(nρ w c w ) ≈ 2.8, thus we expect the water to cool down approximately 3 times as fast as the solid heats up.
The lower row gives the same profile for d p = 5mm and d p = 10mm, and as expected, the heat exchange is significantly slower in these cases.As a consequence, the width of the front increases.
Observe the similarities of the top left and the lower right plots in Fig. 6.This is due to the fact that the thermal transfer coefficient h given by Eq. ( 14) is proportional to d −2 p .For the 2-dimensional radial symmetry case, the two solutions may be proved to be identical up to scalings of t and r.

Energy efficiency
Finally we study the energy recovery from an ATES well based on 12 hours injection and extraction periods.In general, the energy transfer E in the well over a time interval τ is given by where T w,0 = T w (0, t) is the water temperature at the well.The efficiency can me measured in terms of the energy recovery factor given by, (Doughty et al., 1982) Clearly, if the injected water has a constant temperature and the injection rate Q is constant, then E(τ injection ) = ρ w c w (T inj − T 0 )Qτ injection .During pumping the water temperature at the well will vary depending on the dispersion of temperature in the aquifer, which includes natural heterogeneity and the heat transfer between the liquid and the solid phase.Numerical simulations of the temperature of water and solid at the well over five consecutive cycles are given in Fig. 7.The corresponding recovery rates are shown in the lower right corner of the plot.
We observe that the heat exchange has a significant impact on the efficiency rate for d p = 500mm, otherwise not.We notice that the efficiency recovery rate based on this simplified model corresponds well with the rates achieved for an ATES system in the same aquifer presented in Birhanu et al. (2015b).
The temperature profile over one cycle (injection and pumping) is given for the two extreme cases d p = 500mm and d p = 1mm in Fig. 8.

Conclusion
This paper has briefly reviewed the differential equations for heat transport in water-saturated porous media, and presented numerical and analytical solutions for radially symmetric flow.In particular, a simple similarity solution was obtained for the heat transfer in a 2D horizontal confined aquifer in local fluid/solid thermal equilibrium.For a time varying fluid flow and different fluid and solid temperatures, that is, the non-equilibrium or delayed case, solutions have to be obtained numerically.

5
The numerical algorithms have been based on a semi-discrete Lagrangian formulation.
The numerical models have enabled us to consider the primary purpose of this investigation, namely to calculate the recovery factor of a one-well ATES system with a cyclic repetition of injection and pumping.It has turned out that the performance is dependent on the total length of the cycle relative to the time scale for the heat transfer between fluid and solid.The latter may be linked to the typical grain size d p as shown in Table 2.For a total cycle of length 24 hours, referring to Fig. performance is seen to be virtually independent of the grain size as long as d p is less than about 100mm (∆t ≈ 25minutes), but significantly affected for d s = 500mm (∆t ≈ 10hours).In the latter case, the efficiency is also significantly reduced.
Based on the presented results, the analytic and numerical solutions should provide a consistent tool for the understanding of water and solid temperatures near wells with radial flow.
)Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-303Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 8 August 2017 c Author(s) 2017.CC BY 3.0 License.3Analytical SolutionsAnalytical solutions express the relation between the principal variables involved directly.This provides basic insight to the problem without any further numerical evaluation.Besides this intuitive conceptual advantage, classical applications of analytical solutions are practical parameter estimation and sensitivity analysis.In practical hydrology physical parameters like permeability or storage capacity, are not always accessible for direct measurements.Response functions (viz.hydraulic head;

Figure 1 .
Figure1.The 1-D similarity solution Eq. (31) and the exact solution Eq. (32) compared with the hyperbolic front for two different values of the diffusion coefficient λ.

Figure 2 .
Figure 2. The exact similarity solution Eq. (39) compared with the hyperbolic front for two different values of the diffusion coefficient λ in the 2-dimensional radial symmetric case.

Figure 3 .Figure 4 .
Figure 3.The extended domain and the characteristic lines in the injection phase.

Figure 5 .
Figure 5. Temperature distribution of water (left) and solid (right) in the aquifer as a function of time t and radial distance r for v(t) = cos(πt).Hot water is injected for t < 0.5 and pumped for t > 0.5.

Figure 6 .
Figure 6.The temperature profile of the solid and the water for different values of dp[mm] in 2D radial flow near a well.The solid lines indicate the temperature of the water and the broken lines for solid temperature.The upper row emphasize water and solid temperature profile of the same particle size with different timescale while the bottom row emphasize the water and solid temperature profile of different particle size over the same timescale.

Figure 7 .
Figure7.The temperature of the water and solid at the well for five consecutive cycles of 24 hours for aquifers of different particle sizes.In the lower right corner, the corresponding recovery rates.

7 Acknowledgments 5 Figure 8 .
Figure 8. Temperature profiles of water (solid lines) and solid (dashed lines) for different particle sizes The left column shows the temperature profiles during injection, the right during pumping.

Table 1 .
Physical and thermal properties of fluid and aquifer for the thermo-hydraulic modelling of the Gardermoen aquifer.

Table 2 .
The heat transfer coefficient (Eq.14) and the estimated time scale (Eq.16) towards thermal equilibrium for different particle size.