New Lie Symmetries and Exact Solutions of a Mathematical Model Describing Solute Transport in Poroelastic Materials

: A one-dimensional model for fluid and solute transport in poroelastic materials (PEMs) is studied. Although the model was recently derived and some exact solutions, in particular steady-state solutions and their applications, were studied, special cases occurring when some parameters vanish were not analysed earlier. Since the governing equations are nonintegrable in nonstationary cases, the Lie symmetry method and modern tools for solving ODE systems are applied in order to construct time-dependent exact solutions. Depending on parameters arising in the governing equations, several special cases with new Lie symmetries are identified. Some of them have a highly nontrivial structure that cannot be predicted from a physical point of view or using Lie symmetries of other real-world models. Applying the symmetries obtained, multiparameter families of exact solutions are constructed, including those in terms of elementary and special functions (hypergeometric, Whittaker, Bessel and modified Bessel functions). A possible application of the solutions obtained is demonstrated, and it is shown that some exact solutions can describe (at least qualitatively) the solute transport in PEM. The obtained exact solutions can also be used as test problems for estimating the accuracy of approximate analytical and numerical methods for solving relevant boundary value problems.


Introduction
The poroelastic theory considers the system formed by an elastic material with pores that can be penetrated by a fluid and/or dissolved solutes.A typical biological example is the glucose solute penetrating a tissue layer.The basic theory of poroelasticity was developed by Biot [1], and its state-of-the-art can be found in the well-known books [2][3][4] (see also recent studies, e.g., [5][6][7][8][9]).A poroelastic material is considered as the superposition of two continuous media: the matrix (skeleton), occupying the fractional volume θ M , and the system of pores saturated by a fluid, occupying the fractional volume θ F (θ F + θ M = 1).The deformation of the system under the fluid pressure is described by a deformation vector, and the dynamics of the deformation under the forces needs in general to be described by second order tensors.The relationship between stress and strain is usually assumed to be linear.The flux of fluid depends on the hydrostatic pressure and solute gradient (called osmotic pressure).The diffusive and convective transport mechanisms should be taken into account as well.The general three-dimensional theory describing transport in poroelastic materials (PEMs) is very complex because relevant mathematical models involve 3D nonlinear partial differential equations (PDEs).As a result, in order to obtain analytical results, one-dimensional versions are usually discussed [10][11][12].Here, we study a one-dimensional model for fluid and solute transport in poroelastic materials (PEMs) that was derived in [11] and generalised in [12].The governing equations of the model read as where σ 1 = σRT and ρ F = ρ 0 F are positive constants, and the lower subscripts t and x denote differentiation with respect to these variables.The physical/biological meanings of the notations used above are presented in Table 1.
Table 1.Description of the symbols used in Equations ( 1)- (5).The nonlinear system of PDEs (1)- (5) was integrated in the stationary case.As a result, all steady-state solutions were identified, and examples of their application for the glucose fluid transport in a biological tissue were provided [11].In the nonstationary case, system (1)-( 5) is not integrable; therefore, the classical Lie method [13][14][15][16][17] was adopted for search exact solutions.Nowadays, this method is widely used for construction of exact solutions of nonlinear PDEs arising in real-world applications and the most remarkable works are cited in the above cited books.However, it can be easily noted that there are not many studies devoted to multicomponent systems of PDEs because essential technical difficulties occur if one intends to find exact solutions for such systems by applying Lie symmetries.Taking into account the above observation, we refer the reader to recent works [18][19][20][21][22][23], devoted to applications of the classical Lie method to the nonlinear three-component systems of PDEs.

Symbol
Although several nontrivial Lie symmetries were identified and successfully applied for finding exact solutions in [11] (see also generalisations in [12]), some limiting cases have been not analysed therein.In particular, the special cases D = 0 and/or S = 0 were not examined in [11,12].It is proved here that the special cases listed above lead to a rich Lie symmetry involving symmetry operators that do not occur for system (1)- (5) with DS ̸ = 0.It should be stressed that the governing equations with D = 0 and/or S = 0 can still be used as a real-world model in some cases.For example, it is known that the diffusion term in Equation ( 5) is negligible (i.e., D = 0) if molecules of large size (because of high atomic weight) are dissolved.A typical example is albumin, for which the diffusivity is circa 50 times smaller than the glucose diffusivity [24].
This paper is organised as follows.In Section 2, all possible extensions of the Lie algebra of invariance of the nonlinear system of PDEs (1)-( 5) are derived.It is proven that there are four inequivalent cases when the system admits additional Lie symmetries depending on the values of the parameters D and S. Section 3 is devoted to the constructions of exact solutions of system (1)-( 5) with S = 0.Because the system admits an additional Lie symmetry, new reductions to systems of ODEs and, as a result, new exact solutions of the nonlinear model in question are derived.In Section 4, an exact solution is analysed in order to show its applicability for modelling solute transport in PEM.The analysis is supported by 3D plots of the solution.In Section 5, new exact solutions of the nonlinear system (1)-( 5) with arbitrary D and S are constructed.The solutions obtained were not identified in the previous study [11].In Section 6, we present conclusions highlighting the main results obtained.

Lie Symmetries
Here, we start from the governing Equations ( 1)-( 5) of the model for fluid and solute transport in poroelastic materials that were derived in [11].By the application of the substitution where p * is a so-called effective pressure, the above system takes the form In (7), parameters k, λ * , D, and S are constants satisfying the restrictions Theorem 1 ([11]).System (7) with arbitrary given parameters k, λ * , ρ 0 F , D, and S is invariant under an infinity-dimensional Lie algebra generated by the Lie symmetries: where g(t) is an arbitrary smooth function (hereafter, the notations According to the standard terminology, Lie algebra (9) is called the principal algebra of system (7) (see, e.g., Chapter 1 in [17]).Because the latter involves several parameters, a Lie symmetry classification problem arises that was not solved in [11].
We remind the reader that Lie symmetry classification (LSC) problems (group classification problems) are the most difficult of those in Lie symmetry analysis.One may claim that Sophus Lie carried out some theoretical foundations for solving LSC problems; in particular, he has carried out classification of Lie symmetries for a class of linear twodimensional PDEs.During the last few decades, a vast number of papers were published devoted to theoretical foundations and applications of algorithms for classification of Lie symmetries of specific PDEs (systems of PDEs).Unfortunately, a misleading terminology is used in many of them because instead of a complete solving of the LSC problem for a given PDE, only examples of extensions of the principal algebra are presented.The current state of the art and the relevant list of most important publications can be found in [17] (see Chapter 2 therein).Theorem 2. System (7) with restrictions (8) admits the extensions of the principal algebra (9) only in the cases presented in Table 2.

Restrictions
Lie Symmetries Extending Algebra (9) H 1 , H 2 and h are arbitrary smooth functions.
Proof.The proof of this theorem is based on the infinitesimal criteria of invariance, which was formulated by S. Lie in their pioneering works [13,14].In the case of a system of PDEs of arbitrary order, this criteria can be found, e.g., in [15] (see Section 1.2.5).Note, system (7) consists of five PDEs of the second order.So, a linear first-order operator of the form (here, the coefficients ξ 0 , ξ 1 and η i , i = 1, . . ., 5, to-be-determined functions of independent and dependent variables) is a Lie symmetry (operator of Lie's invariance, point symmetry) of system (7) provided that the following equalities are simultaneously satisfied for each solution u(t, x), ρ(t, x), p * (t, x), θ F (t, x), c(t, x) of the PDE system (7).Here, X (2) is the second-order prolongation of the operator X, which is again the first-order operator with coefficients defined by the well-known formulae via the first-and secondorder derivatives of unknown coefficients ξ 0 , ξ 1 and η i (see, e.g., [15], Section 1.2.1).After relevant calculations, Formula (10) are reducible to a linear system of PDEs, called the system of determining Equations (DEs), for finding the functions ξ 0 , ξ 1 and η i .Some of the equations belonging to the system of DEs do not contain parameters S and D; therefore, solving them, we obtain the most general form of the Lie symmetry operator where α i (i = 0, . . ., 3) are arbitrary constants, and g(t) is an arbitrary smooth function, while the constant α 4 and the functions η 4 and η 5 must satisfy the restriction Sα 4 = 0 and the remaining DEs.
If one assumes that all parameters arising in (7) are arbitrary constants, then the result presented in Theorem 1 in a straightforward way is derived.However, here, we need to determine all possible values of parameters (under restrictions (8)) when a wider Lie symmetry occurs.Looking at the system of DEs ( 12)-( 15), one notes (see Equation (12)) that two different cases, S ̸ = 0 and S = 0, should be separately examined.
In the case S = 0, an extension of the principal algebra (9) via the operator Finally, special values S = 0 and D = 0 lead to the widest Lie symmetry.In fact, Equations ( 12)-( 14) simply disappear and the remaining equations are Solving Equation ( 18), one finds Substituting the above function into Equation ( 17), the function Thus, the Lie symmetry operator (11) with the above functions η 4 and η 5 produce exactly the Lie symmetries Although the cases listed in Table 2 are very special and are questionable from the applicability point of view, new Lie symmetries are very interesting.In particular, system (7) with S = 0 admits the symmetry x∂ p * + x 2 2λ * ∂ u that can be not predicted using any physical laws.Here, we study this case in detail and our goal is to construct a wide range of time-dependent solutions.

Ansätze, Reductions and Exact Solutions in the Case S = 0
As it was noted above, the most interesting cases from the symmetry point of view are the second and the fourth cases (see Table 2) because they involve the Lie symmetry 2λ * ∂ u that cannot be predicted by any physical/biological consideration.The fourth case is rather artificial from an applicability point of view.Indeed, the restrictions D = S = 0 means that both diffusive and convective transports are neglected, which is a very strong assumption (however, one cannot claim that such a special case is absolutely unrealistic).Thus, we study the second case, S = 0, in what follows.
The general technique for constructing an exact solution using Lie symmetries is wellknown from the theoretical point of view.There are two different approaches described, e.g., in Section 1.3 [17].One of them is based on the construction of optimal systems of inequivalent (nonconjugated) subalgebras of the given Lie algebra of symmetries.Although this technique is very popular, one successfully works only in the case of Lie algebras of low dimensionality (up to 4).In fact, all optimal systems for such algebras were described in a seminal work [25].However, a pure algebraic problem occurs if one deals with a Lie algebra of the dimensionality 5 and higher.To the best of our knowledge, there is no generalisations of the results presented in [25] on higher-dimensional Lie algebras.Notably, the Lie algebra of symmetries of system ( 7) is infinity-dimensional (see the operator g(t)∂ p * in Theorem 1).
The second technique for constructing exact solution via Lie symmetries is based on the most general linear combination of basic operators of the Lie algebra in question.Of course, many technical difficulties arise in the case of higher-dimensional Lie algebras; however, these difficulties can be overcome by a careful analysis of the relevant invariant surface condition.Here, this is demonstrated in the case of Lie algebra with the basic operators Let us consider the most general linear combination of Lie symmetries listed above where α 1 , α 2 and β i (i = 0, 1, 2, 3) are arbitrary constants, then there is no reduction to an ODE system).
Depending on the value of the constants α 1 and α 2 , operator (19) produces three inequivalent ansätze for the reduction of system (7) to ODE systems.Thus, one needs to consider such cases: From the very beginning, it can be noted that the function g(t) does not play a role in the first two cases, and we may set g(t) = 0 without losing a generality.In fact, ansätze obtained contain the added term g * (t) in expressions for the function p * , and this term vanishes after substituting into system (7).So, one always can generalise an arbitrary solution by inserting into the pressure component the additional term g * (t) that simply follows from the Lie group generated by the operator g(t)∂ p * : In the third case, the function g(t) does play a role.
Substituting ansatz (20) into PDE system (7), one arrives at the ODE system where new notations We are looking for exact solutions of the ODE system (21) in the two essentially different subcases (i) 1 + F = 0 and (ii) Subcase (i).System (21) after the relevant calculations leads to the result where C i are arbitrary constants and constant D > 0.
If β 2 = 0, then the general solution of Equation ( 22) has the form f 5 = C 5 + C 6 ω.Thus, taking into account ansatz (20) and the above formulae, the exact solution of the PDE system (7) with S = 0 (here, θ F is an arbitrary smooth function) is obtained.
If β 2 ̸ = 0, then the general solution of Equation ( 22) can be found in an explicit form only for some fixed function f 4 .Let us solve Equation ( 22) assuming f 4 (ω) = θ 0 ω n (here, θ 0 > 0 and n are arbitrary constants).So, one needs to construct exact solutions of Equation Equation ( 26) has essentially different general solutions depending on the value of the constants β 2 and n.In the simplest case, when n = 0, the general solution of Equation ( 26) has the form The obtained functions f 4 = θ 0 and f 5 from ( 27) lead to the components of the exact solution of the PDE system (7) with S = 0 and D > 0, while other components of the corresponding solution are presented in Formulae ( 23)- (25).Now, let us return to Equation (26).In the case n ̸ = −2, this equation is reduced to the Bessel equation if β 2 < 0, and to the modified Bessel equation Thus, the general solution of Equation ( 26) with n ̸ = −2 can be expressed in the following forms , (28) if β 2 < 0, and Here, J and Y are the Bessel functions, while I and K are the modified Bessel functions.
Note that solution (27) can be derived from ( 28)-( 29) by setting n = 0 in the last ones.
In the case n = −2, the general solution of Equation ( 26) has the form Note that in the case n = −2 and 4β 2 θ 0 + D < 0 Equation (26) has only complex solutions, not real ones.The functions f 5 (from (28)-( 30)) and f 4 = θ 0 ω n produce the components of the solutions of the PDE system (7) with S = 0.The remaining components of the corresponding solutions are still given by Formulae ( 23)- (25).Subcase (ii).In this case, the functions f 1 , f 2 , f 3 and f 4 can be expressed via the function F, namely, where the function F is a solution of the equation while the function f 5 must be founded from the equation (33) ODE (32) is a nonlinear Abel-type equation and its full integration is questionable; moreover, even particular solutions are unknown.However, assuming that the function F is linear, one can find the exact solution provided that the following restrictions are satisfied Another possibility to integrate ODE (32) arises if one sets In this case, ODE (32) takes the form The latter is integrable via the Lambert function, where The exact solution (36) can also be presented in the implicit form (F(ω)) It can be easily seen that a special case occurs when λ * α 2 ρ 0 F = 1.In this case, one again arrives at the linear function Each of the exact solutions (34)-(37) can be used for finding the function f 5 by solving the linear ODE (33).Because the formulae obtained are very cumbersome (in particular, the Heun functions are obtained if one applies (37)), here, we present only the details concerning solution (34).In this case, ODE (33) can be rewritten in the form Although this equation is very awkward, one may find the substitution transforming ODE (38) into the well-known Whittaker Equation (see, e.g., [26]) where The general solution of Equation (40) has the form where M and W are the Whittaker functions.In the case µ = 1 2 + ν (that leads to the conditions β 2 = 0 or β 2 = C 0 α), the function M takes the form M 1 2 +ν,ν (τ) = e − τ 2 τ 1 2 +ν .Thus, the general solution of the Whittaker equation in this special case has the form Taking into account (39) and (42), we arrive at the solution of ODE (38) under restrictions β 2 = 0 and β 2 = C 0 α, namely, Using the obtained functions f 5 and the above formulae, one can construct exact solutions of the PDE system (7).In particular, setting C 6 = 0 in (43) for simplicity, an exact solution of the system in question can be written down in the form where Similarly, setting C 6 = 0 in (44), another solution of the PDE system (7) with S = 0 can be derived in the form while other components have the same form as in (45).
Here, we substitute ansatz (20) with P 2 ̸ = 0 into PDE system (7) with S = 0 to obtain a reduced ODE system.As a result, one arrives at the system The first four equations of the ODE system (47) can be easily integrated, producing the functions while the last equation on the function f 5 (x) takes the form It is difficult to construct exact solutions of Equation ( 49) without additional restrictions.Note that Equation (49) in the case C 4 = 0 is the Heun equation, and its general solution can be constructed using software such as Maple 2017.3.Since this case leads to very cumbersome formulae, we omit those here.
Let us set β 0 = β 1 = β 2 = C 4 = 0, β 3 ̸ = 0, then Equation (49) takes the form The latter is integrable in the terms of modified Bessel functions.Thus, we arrive at the exact solution of the PDE system (7) with S = 0 : Let us set β 0 = β 3 = 0, β 1 ̸ = 0, then Equation (49) takes the form Applying the substitution to Equation (50), we again arrive at the Whittaker Equation ( 40) with the parameters Note that we consider only the special case µ = 1 2 + ν that leads to the condition Now making the same calculations as in subcase (ii), we obtain the general solution of Equation (50) as Setting for simplicity C 8 = 0 and using Formulae ( 46), ( 48) and ( 51), we arrive at the exact solution of the PDE system ( 7) β 1 e β 2 t .
In case I I I, operator (19) produces the ansatz (we remind the reader that the operator g(t)∂ p * should be taken into account in this case): where f i (t) (i = 1, . . ., 5) are new smooth functions.Substituting ansatz (52) into the PDE system (7) with S = 0, we arrive at the condition β 3 = 0.In this case, the corresponding ODE system has the form Solving system (53), we obtain the exact solution of the PDE system (7) with S = 0 as follows where G ′′ = g(t) ρ 0 and p * (t) are arbitrary smooth functions.
Remark 1.In order to find exact solutions, the condition β 3 = 0 was used in several case.It means that the operator x∂ p * + x 2 2λ * ∂ u is not taken into account, hence one may apply the relevant ansätze to system (7) with a nonzero parameter S.This is carried out in Section 5.

An Example of the Exact Solution Describing Solute Transport in PEM
In order to show that a given exact solution describes solute transport in PEM, one should check some properties from the very beginning.In particular, the components ρ(t, x), θ F (t, x), p(t, x) and c(t, x) must be nonnegative; moreover, θ F (t, x) ≤ 1.The displacement function u(t, x) should be either nonnegative (this means that PEM is expanding), or nonpositive (PEM is shrinking).There is also a natural initial condition for the displacement, u(0, x) = 0, i.e., no deformation of PEM in the initial time moment t = 0. Obviously, only some exact solutions satisfy the above requirements on a given space interval [0, L] and for t ≤ T, T > 0.
Here, we look in detail at the exact solution (45).By using the space translation x → x + x 0 (x 0 > 0) and setting , and Obviously, the initial condition u(0, x) = 0 is satisfied.It can easily be shown that all components of the solution (54) are bounded and nonnegative in the domain Ω = {(t, x) ∈ (0, T) × (0, L)} and 0 ≤ θ F (t, x) < 1 if the following restrictions hold : Using the function p * from (54) and Formula (6), one can derive the hydrostatic pressure of the poroelastic materials by the formula Plots of the functions u, θ F , c and p defined in the domain Ω = {(t, x) ∈ (0, 1) × (0, 1)} with the parameter restrictions (55) are presented in Figures 1 and 2. The above restrictions (55) guarantee positive values of displacement, i.e., a given layer of PEM is expanding.
In order to obtain negative values of displacement, i.e., a given layer of PEM is shrinking, one needs the restrictions that are listed below.An example is presented in Figures 3 and 4.
If α > 0, then Thus, we have demonstrated that the exact solution (45) with correctly specified coefficients can describe (at least qualitatively) the solute transport in PEM.

New Exact Solutions of System (7) with S ̸ = 0
In [11], some of the simplest exact solutions of system (7) with an arbitrary parameter S are constructed.Here, we show how new solutions can be derived.As noted above in Remark 1, the condition β 3 = 0 in the ansätze constructed in Section 3 means that those are applicable for (7) with nonzero S as well.
The simplest case occurs in the case of ansatz (46).Substituting the latter into system (7) with an arbitrary parameter S, we arrive at the ODE system Thus, solving the ODE system (57) and taking into account (46) with P 2 = 0, we obtain the solution of system ( 7) , where ρ 0 (x) and θ 0 (x) are arbitrary smooth functions, while g(x) is an arbitrary solution of the linear ODE If β 2 ̸ = 0, then the solutions of the above ODE can be written down in an explicit form provided the function θ 0 (x) is correctly specified (typical examples are (x + x 0 ) k and e kx ). If (here, C 0 and C 1 are arbitrary constants).
Let us construct solutions using restrictions (35).In this case, the function F from (37) can be used to specify the functions f i , i = 1, . . ., 4 in (31).Setting C 5 = −1 for simplicity and using the restriction λ * α 2 ρ 0 F = 1, we arrive at the formulae while the equation for finding the function f 5 takes the form In (58), coefficients A i and B i are defined by the formulae The general solution of Equation ( 58) can be constructed via the Heun functions.To avoid cumbersome formulae, we consider only some cases in which Equation (58) can be reduced to known equations, in particular, the Whittaker equation and the Bessel equation.
In the case of the additional restrictions Equation ( 58) is transformed into the Whittaker Equation (40) with the parameters Taking into account (59), restrictions (60) lead to the condition β 1 = − 2 kρ 0 F .Thus, using Formula (41) and the above substitution, one obtains the function f 5 in an explicit form.Finally, applying ansatz (20) with the above specified beta-s and renaming C 4 → − 2C 4 α 2 kρ 0 F , the following solution of the nonlinear system (7) was constructed: where α ̸ = 0, β 0 , β 2 , p 0 , C 4 and C 5 are arbitrary constants To construct exact solutions of Equation (58) in terms of elementary functions, we consider the special case µ = 1 2 + ν, which leads to the conditions Thus, the general solution of Equation (58) with Thus, the component c(t, x) of the solution (61) of the PDE system (7), corresponding to solution (62) with C 6 = 0, has the form Setting C 6 = 0 in (63), another solution of the PDE system (7) can be derived in the form , while other components have the same form as in (61).
In the case when the additional restrictions Equation ( 58) is transformed into the Bessel equation if B 2 > 0, and to the modified Bessel equation Thus, the general solution of Equation (58) under restrictions (64) takes the form Finally, using ansatz (52) and making similar calculations to those in the end of Section 3, one obtains the following exact solution of the PDE system (7): G ′ (t) .

Conclusions
Here, the one-dimensional model for fluid and solute transport in PEM based on the nonlinear PDE system (7) is studied.Although the model was recently derived and some exact solutions, in particular steady-state solutions and their applications, were studied in [11,12], special cases occurring when some parameters vanish were not analysed therein.Since the governing equations are nonintegrable in the nonstationary case, the Lie symmetry method and modern tools for solving ODE systems are applied in order to construct time-dependent exact solutions.Depending on parameters arising in the governing equations, several special cases with new Lie symmetries are identified.Some of them have a highly nontrivial structure, see Cases 2 and 4 in Table 2, that cannot be predicted from a physical point of view or using Lie symmetries of other real-world models.Applying the symmetries obtained for the reduction in the governing PDEs, multiparameter families of exact solutions are constructed, including those in terms of elementary and special functions (hypergeometric, Whittaker, Bessel and modified Bessel functions).Notably, the computer algebra package Maple was used for solving known ODEs.
A possible application of the solutions obtained is demonstrated in Section 4. It is shown that the exact solution obtained in Section 3 can describe (at least qualitatively) the solute transport in PEM provided that relevant parameters are specified correctly.The obtained exact solutions can also be used as test problems for estimating the accuracy of approximate analytical and numerical methods for solving relevant boundary value problems for the nonlinear PDE system (7).
In Section 5, new exact solutions are constructed for the nonlinear PDE system (7) with arbitrary parameters.These solutions essentially generalize the results derived earlier in [11].
This work can be continued in different directions.For example, a natural generalisation of the model has been developed in [27] by taking into account internal sources/sinks.One may expect that the generalised model also admits new Lie symmetries provided some parameters vanish.Another work could be carried out in order to generalize the model on higher-dimensional cases and study properties of the model derived, in particular, by applying the Lie symmetry method.
Description u deformation vector ρ mass density θ F fractional volume of fluid phase F θ M fractional volume of matrix phase M ρ F mass density of fluid phase F c solute concentration in PEM p mechanical pressure in PEM σ reflection coefficient of PEM R gas constant T temperature λ + 2µ elastic modulus with Lame constants λ and µ k hydraulic conductivity D solute diffusivity in PEM S = 1 − σ sieving coefficient of solute in the PEM