Porous Medium Equation in Graphene Oxide Membrane: Nonlinear Dependence of Permeability on Pressure Gradient Explained

Membrane performance in gas separation is quantified by its selectivity, determined as a ratio of measured gas permeabilities of given gases at fixed pressure difference. In this manuscript a nonlinear dependence of gas permeability on pressure difference observed in the measurements of gas permeability of graphene oxide membrane on a manometric integral permeameter is reported. We show that after reasoned assumptions and simplifications in the mathematical description of the experiment, only static properties of any proposed governing equation can be studied, in order to analyze the permeation rate for different pressure differences. Porous Medium Equation is proposed as a suitable governing equation for the gas permeation, as it manages to predict a nonlinear behavior which is consistent with the measured data. A coefficient responsible for the nonlinearity, the polytropic exponent, is determined to be gas-specific—implications on selectivity are discussed, alongside possible hints to a deeper physical interpretation of its actual value.


Introduction
Membrane gas separation is a technology with a well established market [1] and understanding the flow of gases through membranes is therefore of crucial importance. Researchers are competing in surpassing the Robeson upper bound [2] by developing new materials and fabrication processes [3][4][5][6][7][8][9][10]; the other direction of research aims at a deeper understanding of mechanisms of transport of gases through the membranes.
To properly classify a membrane, two main factors are considered: permeance (gas permeation rate per membrane area per pressure difference) or permeability (permeance per membrane thickness) and selectivity (ratio of permeances or permeabilities for a pair of gases) [11,12]. This approach implicitly assumes at least one of the following: (i) the system is linear in normalized parameters, namely membrane area, membrane thickness, and pressure difference; (ii) any nonlinearities in said parameters will cancel out; (iii) the nonlinearity is negligible in the operating range of the system. This is a perfectly valid approach as long as a false implication is not assumed, that by determining selectivity one is free to consider the measured permeability as a linear function of the pressure difference.
Various approaches of modeling the transport mechanisms inside graphene/GO have been reported. The most common are molecular dynamic simulations with an implementation of pores geometry and energy barriers [33][34][35][36], accompanied by adsorption and direct penetration in terms of statistical thermodynamics [37] and molecular sieving [38]. The pressure dependence is seldom addressed and if so, then only as a linear parameter of the simulation intended to be normalized.
Our measurements of gas permeability of a GO membrane [3] within a manometric integral permeation apparatus [39] show nonlinear dependence of gas permeability on input pressure, invalidating any linear partial differential equation (PDE), or any other modeling technique regarding the pressure as normalizable, as a proper governing equation/model for the gas transport through the membrane.
In this manuscript we aim to (i) report the observation of the above mentioned nonlinearity; (ii) provide a criterion for a PDE/model of gas transport through a GO membrane to be able to describe the nonlinearity; (iii) propose a Porous Medium Equation (PME) as a candidate for a governing equation. PME is rather a class of nonlinear PDEs of a given form-we derive two specific PMEs for pressure and resp. density in Appendix A to provide mathematical context.
We describe the dynamics of the experiment from the point of view of the integral cell of the apparatus and based on law of conservation of mass we reformulate the dynamics in terms of pressure and density inside the membrane. We introduce some reasoned assumptions and simplifications in order to make it possible to study the pressure and density profiles within the membrane in their stationary state, without actually considering any specific PDEs as governing equations.
The profiles are solutions to a boundary value problem (BVP) corresponding to the specific experimental setting. We simulate the observation of nonlinearity by introducing a multiplicative change in input pressure and calculating the impact on the gas permeation rate-only at this moment the two PMEs and their 1-D stationary solutions of said BVPs are introduced. For the sake of completeness we present analytical solutions to BVPs on a finite domain for Dirichlet and Neumann boundary conditions and means of rewriting mixed boundary conditions as either of them.
Evaluating the consequences of the multiplicative change in input pressure results in a convincingly accurate power-law formula capable of describing the experimental data. This is the key argument in favor of considering PME as a valid governing equation for gas transport through GO membranes. We discuss the implications of the nonlinearity on selectivity, as well as a possible interpretation of the polytropic exponent-the coefficient responsible for the nonlinearity of PMEs.

Manometric Integral Permeameter
The main component of the apparatus is a stainless-steel construction of two cells with a membrane in between them, for the 3-D model and technological scheme see Figure 1. The top one (source cell) supplies the gas-it is connected to a buffer tank with a pressure controller. The bottom one (integral cell) is equipped with a capacitive pressure gauge. The membrane is placed on a porous sintered steel plug which acts as a porous support for the membrane. Both halves can be connected via a system of valves to a rotary oil pump so that it is possible to evacuate the device. Both cells are water-tempered at 25°C via total of eight horizontal tubes. The experiment is controlled by six valves labeled V1-V6.
To make a single gas permeability measurement, both cells are evacuated, then a specific gas at specific (controlled) input pressure is supplied into the source cell. As the gas flows through the membrane an increase in pressure inside the integral cell is detected.
After a given criteria is met (time elapsed or gauge maximum limit reached) the experiment is terminated.

Membrane
The preparation method for GO is reported in [3], the specific membrane manufacturing process was described as follows.
A volume corresponding to 82 mg of GO was taken from the stock suspension with a GO concentration of 20 mg/mL. The GO suspension was diluted with distilled water to a volume of 25 mL and sonicated in an ultrasonic tray (37 kHz, 300 W) for 15 min. The diluted GO suspension was placed in a pressurized filter cell equipped with a polyamide filter membrane. Filtration was performed using compressed air at a pressure of 3-5 bar. After filtration, the still wet GO membrane was moved (still with the polyamide support) to a desiccator with a constant 30% relative humidity (provided by saturated MgCl 2 solution at 20°C). In this environment the membrane was stored for at least 24 h in order to balance the moisture in the membrane.

Mathematical Formulation of the Experiment
The mathematical goal is to develop a criterion any PDE must fulfill in order for it to be able to describe the observed nonlinearity. The derivation is done in general setting-PME is not considered as a specific mathematical model for our experiment in this section.
The complete mathematical description of the experiment leads to a (non)linear evolution equation with appropriately chosen initial and boundary conditions. Unfortunately, such problem can hardly be solved analytically. Hence, a few simplifications are adopted (described below including the reasoning) so that analytical solutions are available. The main idea is to split the solution into two parts: (i) solution in time domain in the integral cell (assumed to be linear at the beginning of the experiment) and (ii) solution in space domain in the membrane (assumed to be stationary very soon). The two solutions are then linked via law of conservation of mass.

Linearity of Single Measurement in the Integral Cell
In this section we claim that in the range of detected pressures the measured response can be approximated as a linear function of time.
Measurements as described in Section 2.1 are a standard experimental identification technique of measuring step responses. The available range of input pressure p in is given by the safety limits of gas tubing and ranges from 1 to 4 bar; the range of our pressure gauge is 10 Torr = 1333 Pa. This means that it is possible to measure approximately only the first L = 0.0133, resp. L = 0.003325 of the response. In this range the response is almost a linear function of time, so the only extractable information is its slope-the pressure time difference as a ratio p end t end .
The pressure difference changes slightly during the measurement-to vindicate the negligibility of the effect it can be roughly bounded by comparing a linear function with a standard unit step response of first order system, where the first order system is the only possible long-term-behavior-explaining approximation possible, provided we only the slope is measured. Let the unit step response hit the limit value L, At that time the differences in values and derivatives are which for the greatest achievable limit value L = 0.0133 yields ∆y ≈ 9 × 10 −5 and for the lowest achievable limit value L = 0.0026 yields ∆y ≈ 5.6 × 10 −6 . Furthermore, the reported data never actually got even close to the gauge limit-the maximum value of L is another order of magnitude lower at L = 0.000126, which corresponds to the upper bound of the error being less then than 1‰, which is indeed negligible in comparison with other effects such as the input pressure deviating from the Heaviside function, the precision of membrane thickness measurement, or even the homogeneity of the thickness.

Nearly-Stationary Behavior within the Membrane
In this section we claim that if the measured response is a linear function of time, then the pressure profile inside membrane is close to stationary state.
The measured response of the system should be slightly delayed, as it takes some time for the gas to actually move all the way through the membrane and to the gauge. An initial delay can be observed; however, it cannot be distinguished between a delay caused by membrane, by the gauge distance, and by control timing uncertainties. Rather a few seconds are excluded from the evaluation, until the profile becomes linear.
The fact that the time profile becomes linear and stays linear for the measured time interval hints strongly towards an idea that the inner dynamics of pressure distribution inside the membrane is much faster than the dynamics of filling of the integral cell.
We claim that for sufficiently large time t > ε the shape of the profile in the membrane is nearly static, otherwise a different than linear time dependence would occur; an argument supported also by comparing sheer volumes-the membrane is three orders of magnitude smaller in volume than the integral cell.

Conservation of Mass on the Membrane/Integral Cell Interface
In this section the goal is to express the measured quantity, pressure time derivative inside the integral cell-we do that by comparing the mass accumulation inside the integral cell from the point of view of ideal gas law and the gas flux through membrane.
From the ideal gas law the mass inside the integral cell can be expressed as where the subscript IC stands for property of the integral cell and p IC is the mean pressure. Assuming that only mass and pressure can change, taking the time derivative yields There is only one way of changing the mass inside the integral cell-by gas flowing through the membrane. This provides another way of expressing the mass accumulation, commonly described in terms of surface integral of mass flow through the boundary, The cell is assumed to be perfectly insulated, so the only part of the boundary with nonzero flux is the membrane with surface A M . If we think of the membrane as being perpendicular to the x-axis and localized at x ∈ (0, w), then the integral cell can be localized at x ∈ w, w + V IC A M by introducing a cylindrical simplification of the cell. From the simplification it follows that there is no radial flux anywhere, allowing us to simplify and evaluate the integral, ‹ The limit stresses our point of view of the integral cell still being considered. Based on conservation of mass, we can switch the point of view and state that whatever comes inside the integral cell must have come from the membrane, where all the necessary tools are at hand. Introducing Darcy's law of diffusion, where the left-sidedness of the limit switches the point of view to that of the membrane. Equation (7) can be than rewritten in terms of the properties of the membrane, Combining the two expressions of time change of mass, first derived from ideal gas law (5), second from the surface integral (9) the pressure time derivative can be expressed, The key observation here is that the measured time-dependent increase of pressure is proportional to the density and pressure gradient at the bottom of the membrane.

Boundary Conditions for Pressure Profile in the Membrane
In this section the goal is to interpret the experiment as a boundary value problem, whose stationary solution will be further investigated.
The focus is now on the membrane itself. The boundary condition at x = 0 is fairly simple to implement-the pressure is set and controlled there, thus a pressure Dirichlet boundary condition is prescribed.
As it was already hinted by the near-stationary state argument, we will not be predicting the whole dynamics. We rather prescribe the second boundary condition as static and focus on the consequences of changing the input pressure.
To get the stationary solution an actual condition must be of course prescribed at x = w. To be as close as possible to the actual measurement process, we prescribe a general Dirichlet condition p out with the intention of taking the limit p out → 0 + , if required.

Representation of the Nonlinearity
The final goal is to formalize a method of PDE validation via the stationary solutions' capability of predicting the observed nonlinearity of the relation between input pressure and gas permeation rate.
We will work around the dynamics by exploiting Equation (10)-if two different situations can be described, one for input pressure p in and second for a different input pressure ξ p in , by comparing their values of ρ(w − ) and ∂p ∂x (w − ) the dependence of ∂p IC ∂t on the multiplicative factor ξ can be predicted.
Let F be a function of p in returning the value predicted by (10) for a general p out , where the ρ and ∂p ∂x are derived from the solution to the BVP p(0) = p in , p(w) = p out . There are two possibilities of how to introduce the multiplicative factor-either as a BVP p(0) = ξ p in , p(w) = p out or as a BVP p(0) = ξ p in , p(w) = ξ p out . The former case would be a better representation of an experiment terminated by reaching the gauge limit, the latter case would correspond better to a termination-by-time-elapsed experimental setting. We will solve the latter and at the end show a solution to the former. Comparing the solutions with and without ξ-factor should predict (hopefully nonlinear) dependence, which for linear case would simply read f (ξ) = ξ.

Numerical Simulation Support for the Assumptions
To further support the claims and assumptions, most importantly the nearly-stationary behavior within the membrane, we performed numerical simulations of the membrane and integral cell via PME and found them consistent with both the presence and the shape of the nearly-stationary state within the membrane, as well as with the observed nonlinear dependence of gas permeation rate on input pressure [40].
Using Finite Volume Method the complete dynamics of the system consisting of the membrane and the integral cell was simulated. The geometry was slightly more complicated, as the integral cell was modeled as a coaxial cylinder but with twice the radius of the membrane-the rotational symmetry has been conserved. The single slice being rotated had the inner angle of 2°, five control volumes alongside the membrane radius, and resp. ten control volumes alongside the integral cell radius. Along the common axis the membrane was divided into 50 slides of constant width, the integral cell had 100 slides, linearly increasing their width to fit the actual volume of the integral cell.
The governing equation inside the membrane was derived from the divergence form of general PME, ∂u ∂t = D∆ · Γ|u| Γ−1 ∇u , in the integral cell a linear form of the equation (Γ = 1) with diffusion coefficient of nine orders of magnitude higher has been at place.
For the temporal domain, equidistant Crank-Nicolson implicit scheme with the forward weight of 0.9 and time step ∆t = 0.1 s was utilized. The implicit part of the scheme took only the gradient (∇u) into consideration, as calculating the nonlinear part |u| Γ−1 implicitly would be a disproportionately complex problem to solve.

Stationary Solutions of Porous Medium Equation
Properly derived in the Appendix A, we present three formulations of PME for each quantity in Table 1. PME has been thoroughly studied by Vazquez [41]. In Chapter 4.1-Some very simple solutions he describes 1-D stationary solutions of PME. The actual derivation is rather tedious, highly repetitive, and not of any particular interest, so they are presented in a factual manner.
The mixed BVPs have no explicit solution for general γ, as it contains a root of a fractional polynomial, thus only the system of equations to be solved is presented in Table 2, to be used to rewrite the mixed problem as either Dirichlet or Neumann. Solutions to Dirichlet BVPs are presented in Table 3 and to the Neumann BVPs in Table 4: each solution consists of two profiles, p(x) and ρ(x).
For an in-depth mathematical discussion see [41], particularly Chapter 20.1-Largetime behavior of the GDP, Non-negative solutions for the convergence analysis, and corresponding parts of Chapters 5-8 for the uniqueness analysis. Table 1. Overview of PME formulations. Differential operators act on spatial dimensions, where ∇u is the gradient, ∇ · u is the divergence, and ∆u is the Laplacian.

Form Pressure Density
Laplacian

Nonlinearity Predicted
Once the analytical solutions of stationary BVPs of PME are known, it is possible to calculate the (nonlinear) prediction of the gas permeation rate dependence on input pressure, as described in Section 3.5.
Take the density solution and the derivative of pressure solution to Dirichlet BVP p in and p out , When multiplied together, the whole term with spatial coordinate actually cancels out, making the limit calculation rather trivial, Plugging the evaluated limit into (11), the prediction of the gas permeation rate yields Now compare the prediction after introducing multiplicative factor ξ as in (12), Resorting back to the ambiguity of whether or not to multiply the output pressure by ξ as well, in that case we get which is much harder to cope with. However, as p out is negligible to p in , taking the limit results in exactly the same expression.

General Setting
The volume of integral cell of the apparatus was determined via a combination of geometric calculations and datasheet values as V IC = 42.81 ± 0.005 cm 3 . All presented data was measured on a single membrane which without the polyamide support had a thickness of w = 29.5 ± 0.7 µm and effective area A M of a circle with diameter of 3 cm. Its gas permeability was low enough in terms of pressure gauge limit, thus the criterion for experiment termination was passing of t end = 60 s. In Figure 2 we present normalized measurements of p IC in time and in Figure 3 we present gas permeability [42] dependence on input pressure B(p in ), calculated by Every line, resp. every point is an average of usually four distinct realizations of the experiment.

Measurements and Data Regression
In Figure 2, a slight delay of the system reaction can be observed, due to which the determination of slope has been postponed by 10 s. In case of linear dependence of gas permeability on input pressure all the curves would coincide (due to normalization by p in ); however, there is a clear upward trend. The 4 bar propane measurement is excluded as the membrane started to clog itself. There are two measurements of H 2 , the second is denoted as H * 2 and the latter measurement was performed after the propane clogging. The clogging seems to have changed the internal structure of the membrane in such a way that repeatability has been compromised, thus the measurements of C 3 H 8 and H * 2 should not be compared with the other four. In Figure 3 the nonlinearity of gas permeability dependence on input pressure can be observed directly: after normalization, in the linear case a constant value would be expected, but again an upward trend is present-change in input pressure results in disproportionally greater change of permeability. The dotted line is a nonlinear least square regression of the equation which is the power-law representation of the nonlinearity predicted by PME, as a set-point the gas permeability at p in = 1 bar was used. It results in a pair of parameters-b being the predicted gas permeability at p in and γ being the polytropic exponent.  In Figure 3 the nonlinearity of gas permeability dependence on input pressure can be observed directly: after normalization, in the linear case a constant value would be expected, but again an upward trend is present-change in input pressure results in disproportionally greater change of permeability. The dotted line is a nonlinear least square regression of the equation which is the power-law representation of the nonlinearity predicted by PME, as a set-point the gas permeability at p in = 1 bar was used. It results in a pair of parameters-b being the predicted gas permeability at p in and γ being the polytropic exponent. In Figure 3 the nonlinearity of gas permeability dependence on input pressure can be observed directly: after normalization, in the linear case a constant value would be expected, but again an upward trend is present-change in input pressure results in disproportionally greater change of permeability. The dotted line is a nonlinear least square regression of the equation which is the power-law representation of the nonlinearity predicted by PME, as a set-point the gas permeability at p in = 1 bar was used. It results in a pair of parameters-b being the predicted gas permeability at p in and γ being the polytropic exponent.

Polytropic Exponent Interpretation
For each gas, a pair of parameters (b, γ) is obtained. The polytropic exponent can be further interpreted as a molar heat capacity, resulting in values in Table 5. 6. Discussion

Porous Medium Equation Capability of Describing Observed Nonlinearity
We have presented nonlinear gas permeability measurements in dependence on input pressure. Under some assumptions, most notably the nearly-stationary behavior within the membrane, the observed nonlinearity was described, as can be most clearly seen in Figure 3-any linear PDE would result in a simple constant. In other words, PME turned out to be much more suitable model for gas flow through GO membrane than any linear PDE (such as diffusion equation, for example).
Considering PME as a valid model for gas permeability of GO membranes provides an insight into the inner structure of the membrane, as the membranes should be "porous enough" to be thought of as a porous medium as described by Whitaker in [44].

Nonlinearity Implications for Selectivity
The most significant consequence of the nonlinearity is that for two gases with different polytropic exponent, selectivity as a ratio of gas permeabilities depends on the input pressure used for determining the permeabilities. For selectivity between gases X and Y, and for set-point permeabilities b X , b Y measured at p in , the selectivity dependence can be written as which predicts the selectivity at ξ p in . High selectivity is desired, which corresponds to b X b Y but now the nonlinearity also makes it desirable for γ X γ Y (for ξ > 1). In our data, selectivity tends to decrease with increasing pressure-most noticeably, selectivity based on measurements of H 2 and CH 4 gives α H 2 ,CH 4 (1 bar) = 3.423; however, with increased pressure a significant drop occurs with α H 2 ,CH 4 (4 bar) = 2.417.
By no means are we claiming that our membrane can compete in absolute numbers with state of the art membranes, we mainly want to draw attention to the possible significance of pressure gradient dependence. In industrial applications our results might give insight into the problem of finding an optimal operating point and for researchers the main implication is that points on Robeson plot might not be as static as they seem.

Physical Interpretation of Polytropic Exponent
The molar heat capacities as presented in Table 5 are a direct result of introducing a change of temperature into a system described by the first law of thermodynamics and ideal gas law. Specific interpretation of the molar heat capacities might be possible, if one would describe the membrane from the thermodynamical point of view. Such interpretation could hint to phenomena happening inside the membrane such as sorption; however, further research and measurements are needed.

Pressure Range of the Measurement
Our manometric integral permeameter has an operating range from 1 to 4 bars, possibly to 5 bars if we shrink our safety margins. In this range, only propane measurements exhibited different behavior, clogging the membrane at 4 bars. Other gases might exhibit the same clogging effect at higher pressures-specifying clogging pressures might provide another point of view, e.g., if a correlation with molecular diameter shows up. However, such study would require more experiments and in-depth analysis, which is not in the scope of this paper.

Future Research Outline
We are currently working on reproducing the measured data with a higher resolution of input pressure. Preliminary data show that for a slightly different membrane preparation technique the membranes are at least a bit compressible, which reduces the width but consequently the porosity. To measure the phenomenon we will design a new manometric integral permeameter with a capability of measuring the width of the membrane under pressure to quantify the compressibility and with a higher level of automation.
Introducing compressibility into the equations is the next mathematical step, one which will hopefully result in more insight into the internal structure of the membranes.

Conclusions
Nonlinear dependence of gas permeability on pressure difference of graphene oxide membranes has been reported. A criterion of any model's/equation's capability of describing the nonlinearity has been formulated in terms of a general multiplicative change of input pressure.
Porous Medium Equation has been proposed as a suitable governing equation for the gas flow through a GO membrane based on the analysis of its 1-D stationary solutions on a finite domain and subsequent prediction of the gas permeation rate. The prediction was consistent with the gas permeation rate observed when measuring on the manometric integral permeameter.
Validity of Porous Medium Equation in the measured range provides strong reason to believe that Graphene Oxide membranes are porous in such a way that is required by the definition of the term Porous Medium by [44].
Implications of the nonlinearity were discussed, most notably that selectivity can depend on pressure gradient under which each gas' permeability was measured, if the two gases in concern exhibit different values of polytropic exponent. Data Availability Statement: Processed data is contained within the article, raw measured data is available on demand from L.M. Acknowledgments: Graphene oxide membrane was kindly provided by Dan Bouša from Sofer Group of UCT Prague, the apparatus and gases were kindly provided by Karel Friess, the head of the Laboratory of Membrane Separation Processes of UCT Prague. We are thankful to three anonymous reviewers for their valuable comments and suggestions, which helped to improve the quality of this manuscript.

Conflicts of
posed for a scalar quantity u(x, t), typically in one or three spatial dimensions. He sets the course of derivation for our application; however, he further studies only the general case, as from mathematician's point of view the multiplicative constants can be easily scaled out.
We are more interested in specific forms of PME for pressure and density and the relation between them, thus we derive both of them here, providing useful insight.

Appendix A.1. Prerequisites
To set some ground for the derivation we state the continuity equation and Darcy's law of diffusion. The third equation, which due to lack of nomenclature consensus we will call a polytropic state equation, we derive properly as it is the source of nonlinearity in PME and provides a possible physical interpretation. These three together form the foundation of the derivation of PME.
Continuity equation in porous media, especially the introduction and formalization of medium porosity as the mean ratio of void and solid regions in a control volume, has been thoroughly established by S. Whitaker [44]. We will use the formulation where ρ is the density of gas, v is the velocity vector field, ε is the porosity of membrane, and ∇· is the divergence operator. Darcy's law of diffusion describes the flow of fluid through porous medium. It has been theoretically derived from Stokes equation [45], for consistency we once again refer to the work of Whitaker [46]. We will use the formulation where µ stands for the dynamic viscosity of gas, k is Darcy's permeability, and p is the gas pressure. Polytropic state equation follows from the first law of thermodynamics and state equation of ideal gas, accompanied by Mayer's relation and exploitation of some properties of ideal gas.
We begin with the first law of thermodynamics in differential form, where δQ is the heat obtained by the system, dU is the inner energy differential, and δW is the work done by the system on its surroundings (we follow Clausius sign convention). We introduce a change of temperature dT in our system. The change of inner energy of ideal gas is dU = C V dT, where C V is the heat capacity at constant volume. The total heat given to the system introduces some unknown heat capacity C, δQ = CdT.
For the work term only pressure-volume work is considered, Substituting back to (A4) we can express the temperature differential, Let us do the same for the state equation of ideal gas. For R being the universal gas constant, in differential form we have Combining (A8) and (A9) in a single equation yields Firstly we scale out n by introducing molar heat capacities c = C n and c V = C V n , pdV n(c − c V ) = Vdp + pdV nR .
Now recall Mayer's relation: c p = c V + R for c p being the molar heat capacity of isobaric process. We can now rewrite the differential form as a separable differential equation, which can be solved for p, At last we introduce density via V = m ρ , where m, as a conserved quantity, is constant. The term Km Porous medium equation (PME) is derived from these three laws: continuity Equation (A2), Darcy's Law (A3), and polytropic state Equation (A15). We will derive PME in terms of pressure and density in three different forms: • Laplacian form, where the Laplace operator acts on nonlinear expression of our quantity-in this form we can easily compare PME to linear heat equation; • Divergence form, which is numerically treatable with Finite Volume Method; • Linearized form, where both differential operators (Laplacian and gradient) act on linear expression of our quantity-in this form the equation can be treated with Finite Differences Method. and the polytropic state Equation (A15) remains as it is, p = e 0 ρ γ . (A31)