Surface Permeability of Particulate Porous Media

The dispersion process in particulate porous media at low saturation levels takes place over the surface elements of constituent particles and, as we have found previously by comparison with experiments, can be accurately described by superfast nonlinear diffusion partial differential equations. To enhance the predictive power of the mathematical model in practical applications, one requires the knowledge of the effective surface permeability of the particle-in-contact ensemble, which can be directly related with the macroscopic permeability of the particulate media. We have shown previously that permeability of a single particulate element can be accurately determined through the solution of the Laplace–Beltrami Dirichlet boundary value problem. Here, we demonstrate how that methodology can be applied to study permeability of a randomly packed ensemble of interconnected particles. Using surface finite element techniques, we examine numerical solutions to the Laplace–Beltrami problem set in the multiply-connected domains of interconnected particles. We are able to directly estimate tortuosity effects of the surface flows in the particle ensemble setting.


I. INTRODUCTION
Liquid transport in particulate porous media, such as sand, is customarily classified into fully saturated, funicular and pendular regimes of spreading [1][2][3][4] . The first two regimes of the liquid dispersion occur at relatively high saturation levels s > s c ≈ 10%, where saturation s is defined as the ratio of the liquid volume V L to the volume of available voids V E in a sample volume element V , s = V L V E . At high saturation levels, above the critical value s c , liquid transport takes place in the pore space either fully or partially filled by the liquid.
Our prime concern here is the special case of liquid dispersion at low saturation levels. As the saturation level drops below the critical value, s ≤ s c , that is to the value relevant to the pendular regime of spreading, the liquid volumes in the porous matrix become isolated [2][3][4] .
As a result, at low saturation levels, the liquid is only contained in the pendular rings formed at the locations of the particle contacts and on the particle rough surfaces, and the liquid transport can only occur over the matrix surface elements, as is illustrated in Fig. 1.
The analysis of this regime of wetting, which is crucial for studies of biological processes and spreading of non-volatile liquids in arid natural environments and industrial installations, has shown that the liquid dispersion has many distinctive features and can be accurately described by the so-called superfast non-linear diffusion equation 5,6 .
Our main concern here is the wetting cycle, when the liquid spreads over a dry porous matrix or over a matrix with a very low background saturation level up to s r ≈ 2%. These conditions are similar to those in the case studied previously experimentally and theoretically in 6 . The main driving force of the dispersion process, as is often the case during the wetting cycle, is capillary pressure developed at the moving front in the process of wetting of dry porous matrix, while the liquid bridges play a role of variable liquid reservoirs of uniform surface curvature. Note that, unlike at high saturation values s > s c , at low saturation s < s c the capillary pressure is developed on the length scale of the surface roughness.
Theoretically, the superfast non-linear diffusion equation belongs to a special class of mathematical models. Unlike in the standard porous medium equation 7 , where the coefficient of diffusion D(s) is either constant (leading to normal, Gaussian dispersive behaviour) or a monotonically vanishing function of saturation (D(s) ∝ s m , m > 0, leading to hypodispersive behaviour), in this special case, the non-linear coefficient of diffusion demonstrates divergent behaviour as a function of saturation, D(s) ∝ (s − s 0 ) −3/2 , leading to hyper-2 dispersive character of the spreading process. Here, s 0 is some minimal saturation level (s 0 ≈ 0.5%), which could be only achieved in a state when the liquid bridges cease to exist completely [3][4][5][6] . Note, in that respect, that in the domain of spreading liquid bridges are supposed to never vanish, so that the condition s > s 0 is always fulfilled in the model, and there is no actual singularity of the mathematical description 5,6 . The character of the dispersion process may often change from hypodispersive to hyperdispersive depending on the spreading scenario, that is on the driving forces of the spreading process [8][9][10][11][12] , especially in multi-phase systems during forced imbibition cycles used in petroleum engineering 12 . If the forcing, external pressure at the system inlet is strong enough, the low saturation zone may not have a chance to be developed, and hypodispersive character of the spreading front is observed. Returning to our mathematical model at low saturation values, in the macroscopic approximation, that is after averaging over some volume element containing many particles of the porous medium, the diffusion process in the slow creeping flow conditions can be 3 described by the following non-linear diffusion equation

Grain
where The details of derivation of (1) can be found in 5,6 , here we note that, the resultant governing non-linear equation (1) directly follows from the conservation of mass principle where φ is porosity defined as φ = V E V , which is further assumed to be constant, and Q is the macroscopic flux density. The macroscopic flux density Q is defined in such a way that the total flux through the surface of a macroscopic sample volume element is given by the surface integral Q · n dS, where n is the normal vector to the surface of the sample volume element.
To obtain (1) from (2), one needs to apply the capillary pressure-saturation relationship 5,6,13 dictated by the liquid bridges behaviour and the local Darcy's law 14,15 describing the surface flow in the rough layer of the particle Here, Nc π , N c is the coordination number, that is the average number of bridges per a particle, p 0 = 2γ R cos φ c , γ is the coefficient of the surface tension of the liquid, φ c is the contact angle made by the free surface of the liquid bridge with the rough solid surface of the constituent particles, R is an average radius of the porous medium particles, q and u are the averaged local flux density and pressure in the rough surface layer, µ is liquid viscosity and κ m is the local coefficient of permeability of the rough surface, which is proportional to the average amplitude of the surface roughness δ R , that is the width of the surface layer conducting the liquid flux 4 One needs to emphasise here that two levels of averaging are involved in obtaining the final governing equation (1). While equations (1), (2) and (3) are 'truly' macroscopic, that is obtained by averaging using a volume element V containing many grain particles, equation (4) is only an average over some rough area of a single particle containing many surface irregularities, so that quantities q and u are also only local averages over that sample surface area.
Therefore, to transit from (4) to the macroscopic description, the spatial averaging theorem formulated in 16 should be applied. That is, using intrinsic liquid averaging . . . l = where V l is liquid volume within the sample volume V , one has u l = p and q l Se S = Q. Here, S is the surface area of the sample volume V with the effective area of entrances and exits S e . Note, the ratio S e /S is not just a geometric property, but also takes into account the connectivity of the porous elements. For example, the effective area of entrances and exits S e is only defined by the pathways open to the flow.
As a result of the two-level averaging where K(s) = κ m Se S is the coefficient of permeability defined by The global surface permeability of the particles K as a function of saturation is one of the main elements of the model to accurately represent liquid dispersion at low saturation levels. It is fully defined by the particle geometry and the geometry of the liquid bridge contact areas, Fig. 1 and Fig. 2.
In particular, the disposition and the size of the liquid bridges on the particle surface, that is the size of the domains Ω 1,2 and the angle α, should play a leading role in defining the resistance to the surface flow. It is not difficult to discern that any variations of the contact area covered by the liquid bridges (pendular rings), that is areas Ω 1,2 shown in Fig.   2, or the value of the bridge volume, should affect the global permeability.
Previously, we have shown that permeability of a single particle element can be determined by means of a solution to the equivalent Laplace-Beltrami boundary value problem formulated in the flow domain Ω 0 with the boundaries ∂Γ 1,2 in Fig. 2 17 . We briefly formulate that problem and summarise the previous results in the next part. Here we note 5 that, based on the analysis of the problem, we have been able to show that in a special azimuthally symmetric case of spherical particles, when the two areas covered by the liquid bridges, domains Ω 1 and Ω 2 in Fig. 2, are oriented symmetrically to each other, that is at α = π, the permeability K is supposed to follow the scaling We have studied several generalisations of the symmetric problem, such as arbitrary oriented domains, α = π, on the surface of the spherical particle, and particles of arbitrary shapes emulating the shape of a real sand grain. While variations of the particle shape was found to produce a relatively modest effect on the particle surface permeability, the orientation of the boundaries, emulating tortuosity effects, was found to produce a stronger impact due to the substantial variation of the distance, on average, between the boundary contours ∂Γ 1,2 . It became clear that while the previously obtained scaling was a good first step to estimate the surface permeability of particulate porous media, a more general case of an ensemble of interconnected particles should be analysed to enhance the model predictive power and at the same time to estimate rigorously the effects of tortuosity of the surface flow in the particle assembly. In this study, we will simulate a general case of an ensemble of many particles linked by liquid bridges. We will concentrate on the bunch of spherical particles, but of different radii and randomly arranged in configurations. We compare the random pack configuration results with some symmetric case to estimate the effects of tortuosity and formulate practical recipes to apply the super-fast diffusion model.

II. MICROSCOPIC MODEL OF THE SURFACE PERMEABILITY OF THE ELEMENTS
Microscopically, the liquid creeping flow through the surface roughness of each particle can be described by a local Darcy-like relationship (4) between the surface flux density q and averaged (over some area containing many surface irregularities) pressure in the grooves u 14,15 . Assuming incompressibility of the liquid and that the liquid layer thickness is constant Equation (4) taking into account (6) can then be transformed into the Laplace-Beltrami 6 equation defined on the surface Γ of the particle Here, ∆ Γ designates the Laplace-Beltrami operator, which is defined on the surface element Γ through the surface gradient ∇ Γ tangential to the surface. Formally, let n Γ denote the unit normal to the surface Γ, Fig. 2. Then, one can define the surface gradient of a smooth function u as ∇ Γ u := ∇u − (∇u · n Γ )n Γ and then the Laplace-Beltrami operator is defined The second assumption δ R = const implies that the surface layer is fully saturated, that is its content is not changing on the particle surface. The approximation of the fully saturated rough surface layer is well fulfilled, if the characteristic pressure amplitude |u| is less than the capillary pressure amplitude defined on the length scale of the surface roughness δ R , which is of the order of δ R ∼ 1 µm in typical sands 18 , as is demonstrated in 14 . That is, |u| < u c = γ δ R , and, for example for water (γ = 72 mN/m) at δ R = 1 µm, this results in |u| < 7.2 × 10 4 Pa.
Alternatively, if the surface layer somehow is not fully saturated, parameter δ R should be interpreted as the characteristic width of the liquid layer within the rough surface layer and one needs to presume that variations of the pressure |δu| are negligible |δu| u c . This is usually the case in slow, creeping flow conditions in porous media, and in fact, it is a criterion for the use of macroscopic approximation to such flows 1 . As is shown in 6 , strong negative capillary pressure on the level of u c are only expected at the moving front, so that the approximation is well fulfilled in the macroscopic flow domain. Note also that, it is always assumed throughout this study that that is the amplitude of the surface roughness (or the width of the liquid layer) is always much smaller than the particle size.
A. Permeability of a single particle element Consider, as the simplest example, a spherical particle of radius R with a closed surface Γ, which is split into three sub-domains Ω 0 , Ω 1 and Ω 2 with the surface boundaries between them ∂Γ 1 and ∂Γ 2 , as is shown in Fig. 2. The location of the sub-domains Ω 1 and Ω 2 to each other on the surface is fixed by the tilt angle α. The sub-domains Ω 1 and Ω 2 correspond to the contact area covered by the liquid in the bridges, while the surface flow, described by (4), takes place in Ω 0 . Our prime concern is permeability of the surface elements, so that we only consider steady state problems.
The distribution of liquid pressure u, as it follows from (7), should satisfy the Laplace- Beltrami equation now defined on the surface of the sub-domain Ω 0 Note that, in fact, the condition of the fully saturated surface layer is not essential in calculation of the flows over one particle element of the porous media. It is sufficient to presume that the variation of the capillary pressure on the length scale of the particle |δu| is negligible, that is |δu| u c . In the case when the surface layer is not fully saturated, parameter δ R should be interpreted as the effective thickness of the layer filled by the liquid. At the same time, liquid pressure variation in the bridges is negligible in slow, creeping flows in comparison to that in Ω 0 . So that, one can assume that 8 which are the boundary conditions to the Laplace-Beltrami Dirichlet boundary value problem. The Dirichlet boundary value problem (8)-(9) has at least a unique weak solution, if the domain Ω 0 and the boundaries ∂Γ 1,2 are smooth [24][25][26] , which, if it is found, allows to calculate the total flux through the particle element where n s is the normal vector to the domain boundaries ∂Γ 1,2 on the surface, δ R is the average amplitude of the surface roughness, that is the width of the surface layer conducting the liquid flux and the line integral is taken along a closed curve in Ω 0 , for example the boundary ∂Γ 1 .
If the total flux Q T is determined, one can define the global permeability coefficient of a single particle K 1 . This can be done, if we assume that the particle has a characteristic size D and so that it can be enclosed in a volume element V = D 3 with the characteristic side surface area D 2 . Then, the effective flux density Q can be represented in terms of K 1 (and the total flux Q T ) if the flow is driven by the constant pressure difference U 2 − U 1 applied to the sides of the volume element.
B. Surface permeability of a sphere in the case of azimuthally symmetric domain boundaries Consider now a spherical particle in an azimuthally symmetric case, when the domain boundaries ∂Γ 1 and ∂Γ 2 are oriented at the reflex angle α = π and have a circular shape. We use a spherical coordinate system with its origin at the particle centre and the polar angle θ counted from the axis of symmetry passing through the centre of the circular contour ∂Γ 1 .
In this case, the Dirichlet boundary value problem (8)-(9) admits an analytical solution, so that particle permeability can be determined explicitly. Indeed, problem (8)-(9), if we assume that the liquid pressure distribution u is a function of θ only and independent of the azimuthal angle, is equivalent to 9 with the boundary conditions The analytic solution to problem (12)-(13) after applying the boundary conditions can be represented in the following form where One can now calculate the total flux and the permeability, using its definition (11), So that, taking D = 2R, Parametrically, the coefficient of permeability (16) is inversely proportional to the particle radius R, so that larger particles create stronger resistance to the flow. Noticeably, the coefficient demonstrates strong dependence on the surface layer thickness δ R , that is K 1 ∝ δ 3 R since it is anticipated that κ m ∝ δ 2 R , so that evaluation of this parameter in applications is crucial for the accurate estimates of the liquid dispersion rates.
One can see, if we take θ 1 = θ 0 , in fact assuming small variations of the bridge size and the pressure over one particle diameter, and θ 0 1, in fact considering small values of saturation, s 1, that the permeability coefficient K 1 tends to zero, that is How does the result affect the super-fast diffusion model (1), and basically how can it be incorporated into the main diffusion equation? If we approximate the permeability coefficient K by K 1 obtained in the azimuthally symmetric case at θ 1 = θ 0 , (17), and, using an approximate relationship between the radius of curvature R sin θ 0 of the boundary contour ∂Γ 1 and the pendular ring volume 2 , one can show that Therefore, finally As it follows from (18), the distinctive particle shape results in logarithmic correction to the Apparently, the correction will mitigate to some extent the divergent nature of the dispersion at the very small saturation levels s ≈ s 0 , smoothing out the characteristic dispersion curves.

C. Surface permeability of a chain of spheres in the case of azimuthally symmetric domain boundaries
Consider now how the problem can be formulated in the case of several particles arranged in a single chain, as is illustrated in Fig. 3 in the case of two coupled by the bridge particles.
To create the flow in the system of two coupled particles, one can set pressure difference between ∂Γ 2 are 'internal', that is common to the bridge linking the flow between the two particles. Apparently, the pressure is supposed to be the same on the two contours and due to conservation of mass in steady state conditions in the absence of sinks and sources of the liquid one has where n s 1 and n s 2 are the outward tangential normal vectors to the boundary contours ∂Γ (2,1) 1,2 , and u 1 and u 2 designate distribution of pressure on each particle respectively. As a result, the problem to define the flow and the permeability of the system corresponds to a system of two Laplace-Beltrami equations and 1 sin θ but with a slightly different set of the boundary conditions and where θ 0 , as before, defines the size of the bridge footprint on the particle surface in the spherical coordinate system with the axis of symmetry passing through the centre of the bridge area, Fig. 3. Since we assumed, due to relatively small variations of pressure over a few grain particles, that all bridges are roughly identical, we have only one parameter θ 0 to describe the bridge size.
Apparently, equations (21) and (22) can be integrated twice, similar to the previous problem of a single particle (12), to obtain where C 1,2 and B 1,2 are free constant parameters to be found from the boundary conditions.
It is not difficult to see from (26), that one has C 0 = B 0 implying continuity of the contact flux. Applying the remaining boundary conditions (23)-(25), from (27) and (28) where Ψ One can now calculate total flux and define permeability of the coupled spherical particles where D is the characteristic length scale of the cross-section in the problem, κ m is local permeability of the surface layer, δ R is the layer width and µ is liquid viscosity.
So that, taking simply D = 2R, One can see that the permeability of a system of two coupled particles K 2 is identical to that of a single particle (16), basically from (16) and (32) 13 It is not difficult to discern by deduction that in a general case of N coupled particles in a chain Note, a setup of many beads coupled by liquid bridges can be potentially used in microfluidics to create flexible water channels 19 . If the radius of curvature of the particle chain is much larger than the particle size, the transport through such a microfluidic system should be defined by the permeability of a single particle, relationship (16), if the particle shape can be approximated by a sphere. One can conclude in this part, that if the porous media configuration is made of parallel chains of particles oriented symmetrically to each other, and the flow is generated along the chains, the surface permeability given by (16) is the exact result.

III. SURFACE PERMEABILITY OF A RANDOMLY PACKED PARTICLE ENSEMBLE
In real systems, the particles are interconnected randomly, so that the effects of tortuosity should substantially affect the permeability of the system 1,20-22 . To analyse those effects, we consider an ensemble of spherical particles randomly packed, as is shown in Fig.   4. The randomly packed configuration of approximately 3000 − 7000 particles has been generated by means of a molecular dynamics technique by applying a constant force to every particle placed in a box with reflecting boundaries (in the perpendicular direction to the box side), and interacting via the Lennard-Jones potential with different characteristic length scales R distributed normally, that is with the probability of the particle radius at ∆R/R 0 = 0.3. In this study, there were particles with three different characteristic dimensions R 1 = 1.3 R 0 , R 2 = R 0 and R 3 = 0.7 R 0 . The resultant porosity in the configurations was about 48%.
To obtain the configuration, the particle temperature controlled by the thermostat has been gradually reduced to bring the system to a minimum energy, frozen state. A representative sample volume element with dimensions L B x , L B y , L B z then was cut off the system, as is illustrated in Fig. 4, containing N S = 13 − 17 particles, see Table I  The Laplace-Beltrami method then has been applied after establishing the position of the liquid bridges coupling the particles in the sample. Two particles (of radii R 1 and R 2 ) are assumed to be coupled by a liquid bridge if the distance between their centres r was only slightly larger than the sum of their radii The size of a single liquid bridge footprint H B on the particle surface can be characterised, as before, by the polar angle θ 0 in the polar coordinate system with the symmetry axis passing through the centre of the circular contour, the boundary of the area covered by the bridge, as is shown in the symmetric case in Fig. 3. That is, H B = 2R k sin θ (k) 0 . Due to the specific geometric properties of the pendular rings (constant mean curvature surface), we Ω 0 Indeed, when s 1, the pressure in the pendular ring p is defined by the smallest radius of curvature r 1 , Fig. 6, p ≈ −γ cos φ c /r 1 , which is related with the second radius so that when s 1, one has r 2 r 1 . Obviously, r 2 defines the size of the area covered by the bridge, H B = 2r 2 = 2R k sin θ (k) 0 . If we have two particles of different radii, say R 1 and R 2 , in contact, the size of the bridge area will be approximately the same r at low saturation levels, s 1, with the difference being proportional to r 1 , that is tessellation of the domains, as is shown in Fig. 7 for one particulate element with two boundary contours. The numerical method has been validated against analytical solutions previously demonstrating prescribed order of accuracy and numerical convergence, see details in 17 .
The number of particles in the sample volume element was negotiated between computational efficiency of the surface finite element method (so that, practically, any mesh resolution can be used to deal with any details on the boundary contours ∂Γ As one can see, problem (21) - (26) and hence total flux Q T through a particle or a chain of particles, (15) or (31), are invariant under the transformation of the particle dimension R provided that the angular size of the bridge θ 0 is fixed. In what follows, we change to non-dimensional description by normalizing length scales by the average radius R 0 of the particles in the sample and pressure by the characteristic capillary pressure p 0 = 2γ cos φ c /R 0 . The flux Q T will be normalised by the characteristic value inspired by the analytical result (15) and by the non-dimensional sample box surface area x,y,z /R 0 andŪ 1,2 = U 1,2 /p 0 . The latter normalization allows to bring simulation results in slightly different geometric settings, as is detailed in Table I, into equivalent conditions suitable for comparison, that is basically     tion should not render any excessive (in excess of the statistical errors due to the particle number fluctuations) influence upon the results, that is the value of the total flux and the 'macroscopic' permeability. A posteriori, one can see that this seemed to be the case, Fig.   8, as in different configurations, Table I, the resultant curves are close and parallel to each other.

Parameters of the configurations
There are two main questions, we would like to answer in this part of the study. First, how does permeability of the particle sample depend on the composition? Basically, how strong are there fluctuations? Secondly, what is the contribution of the tortuosity effects? To obtain statistically meaningful results, we consider several randomly generated configurations, as is summarised in Table I. We would like to stress here, that all configurations have been cut off from statistically independent particle distributions generated with the help of random initial distributions of larger number of particles, as we have described.
As before, we are going to find a weak solution to a system of the Laplace-Beltrami equations ∆ Ω (k) 0 u k = 0 defined on each particle domain Ω we set up continuity conditions, for example on ∂Γ 3 and ∂Γ While on a few external boundaries, Dirichlet boundary conditions are set.
The numerical solution allows to calculate the total flux through the system by summing up the fluxes passing through the external contours, where the Dirichlet boundary conditions are set, either at the top of the pack or at the bottom using (10). The results are summarized in Fig. 8.

20
Remarkably, the reduced flux Q T /S 0 Q 0 as a function of where H B is the bridge size H B = 2R 0 sin θ 0 , behaves linearly in all configurations. This behaviour mirrors the flux dependence observed in azimuthally symmetric analytical solutions, see (15) or (31). The variations in the dependencies between different configurations are observed to be well within the statistical error expected in this case, error bar in Fig 8. At the same time, a comparison with a similar, but a regular arrangement, as in Fig. 3 at R = R 0 demonstrates that there is a clear cut contribution from the effects of tortuosity, solid line in Fig. 8.
Indeed, given identical porosity (φ ≈ 50%) and mean particle size (R/R 0 = 1) in the regular, symmetric and randomly generated configurations, the normalised flux values differ by a factor of two τ s ≈ 2, which can be regarded as the tortuosity of the surface flow.
The value of the surface flow tortuosity τ s ≈ 2 is consistent with the tortuosity values obtained in porous media in different conditions and configurations 22 . For example, both hydraulic τ h and diffusive τ d tortuosity estimated in unsaturated porous media using different permeability models (often used in applications, for example 27,28 ) was found in between 1.5 ≤ τ h,d ≤ 2 at φ = 50% 22 .
It is important that the result, that is the ratio of the total flux in the random and regular configurations does not practically depend on the size of the contour H B , basically the size of the liquid bridge, and hence the value of saturation in the porous media. This implies, that the observed effect is purely down to the distribution of contacts between the particles, but not the particular pathway on each single particle surface. That is, fundamentally, tortuosity in the surface diffusion processes is a geometric factor independent of the particular surface flow regime. At the same time, the pathways, on average, of course, does depend on the bridge size value H B leading to smaller permeability as the size of the contact area diminishes.
This trend is expected, but essentially, the correction to the effective coefficient of diffusion D(s) ∝ K(s) (s − s 0 ) 3/2 , K(s) = τ s K 1 (s) is only down to a single universal factor of two, τ s , representing the tortuosity effects in surface diffusion in particular porous media at low values of saturation, where K 1 (s) ∝ | ln(s − s 0 )| −1 defined by (18) is the surface permeability of a single porous matrix element.
Note, the value of τ s ≈ 2 is also in agreement with experimental observations and a comparison of the super-fast diffusion model with the data, where the tortuosity effects were estimated to reduce the effective permeability twofold 6 .
This is the main result of this study, which can be used in practical applications to calculate permeability in particular porous media. Basically, as the first step, one can calculate permeability of a single, representative element of the media or several elements to obtain some mean value and its dispersion. This way, permeability K(s), via (18), and the diffusion coefficient D(s) for the macroscopic model can be established in the first approximation. Macroscopic permeability K(s) or the diffusion coefficient D(s) then should be corrected by the universal factor of two, τ s , in the macroscopic diffusion model.

IV. CONCLUSIONS
We have demonstrated that the Laplace-Beltrami method can be used to obtain permeability of particulate porous media at low saturation levels and to estimate contribution from the effects of tortuosity. Essentially, analytical results obtained using azimuthally oriented coupled particles can be used with a universal correcting prefactor τ s , as in (39), to estimate permeability of particle ensembles. That is, from the practical point of view, results obtained by analysing single representative element of particulate porous media can be translated into permeability of a particle composition.
PS was supported through the Royal Thai Government scholarship. TP was partially supported through the EPSRC grant EP/P000835/1.

AUTHOR CONTRIBUTION STATEMENT
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.