Spontaneous imbibition in fractal tortuous micro-nano pores considering dynamic contact angle and slip effect: phase portrait analysis and analytical solutions

Shales have abundant micro-nano pores. Meanwhile, a considerable amount of fracturing liquid is imbibed spontaneously in the hydraulic fracturing process. The spontaneous imbibition in tortuous micro-nano pores is special to shale, and dynamic contact angle and slippage are two important characteristics. In this work, we mainly investigate spontaneous imbibition considering dynamic contact angle and slip effect in fractal tortuous capillaries. We introduce phase portrait analysis to analyse the dynamic state and stability of imbibition. Moreover, analytical solutions to the imbibition equation are derived under special situations, and the solutions are verified by published data. Finally, we discuss the influences of slip length, dynamic contact angle and gravity on spontaneous imbibition. The analysis shows that phase portrait is an ideal tool for analysing spontaneous imbibition because it can evaluate the process without solving the complex governing ordinary differential equations. Moreover, dynamic contact angle and slip effect play an important role in fluid imbibition in fractal tortuous capillaries. Neglecting slip effect in micro-nano pores apparently underestimates imbibition capability, and ignoring variations in contact angle causes inaccuracy in predicting imbibition speed at the initial stage of the process. Finally, gravity is one of the factors that control the stabilisation of the imbibition process.

SCiEntifiC RepoRTS | (2018) 8:3919 | DOI: 10.1038/s41598-018-21002-y Whitesides 19 studied capillary rise in noncircular capillaries. Yu 20,21 introduced the fractal method to describe the tortuousness of capillary. Cai [22][23][24] introduced fractal geometry to modify the LW equation in a single tortuous capillary tube and porous media and found that the rise time exponent is related to tortuosity fractal dimension. Li et al. 25,26 studied meniscus rise in branch-like capillaries and tree-like network systems, whereas Shen et al. 27 further applied this concept in shale analysis. Meng et al. 28 studied mechanisms of spontaneous imbibition and the impact of boundary condition, fluid viscosity and wettability. Shou et al. 29 studied geometry-induced asymmetric capillary flowing behaviours, such as multi-section porous layers and trapezoidal porous media.
In micro-nano tubes, dynamic contact angle and slip conditions are two critical characteristics. Martic et al. 30 showed that the contact angle in the LW equation should be dynamic. Contact angle is related to flow velocity [31][32][33] . Several dynamic contact angle models have been derived. One of the most common is the power law model 34,35 , which is recommended and used to simulate experimental data 36 . Another one is the power series model, which was derived by Blake et al. 37 . Hilpert [38][39][40][41][42] derived (semi)-analytical solutions to liquid infiltration into horizontal, infiltration and downward tubes which considers the dynamic contact angle based on the power law and series models. Meakin et al. 43 summarized multiphase fluid flow and reactive transport in fractured and porous media including complete dynamic behaviour of contact lines and contact angles. Petrov et al. 44 formulated qualitative criteria for applicability of hydrodynamic, molecular-kinetic and molecular-hydrodynamic theories on the basis of the wetting-dewetting asymmetry of the dynamic contact angle. Slippage in micronanopores is another important characteristic. Gas slip is obvious and well-studied in both mechanism and application in gas production [45][46][47][48] . Liquid also slip in nanopores. Several authors [49][50][51] have shown that liquid infiltration in micronanopores is significantly different from the ordinary no-slip boundary. Commonly, non-wetting surfaces tend to cause non-wetting liquid slip 52,53 . Wu et al. 54,55 established a universal model for water flow through nanopores, flow capacity of confined water is 10 −1~1 0 7 times that calculated by the no-slip Hagen-Poiseuille equation for nanopores with various contact angles and dimensions 54 . Sometimes wetting surfaces also have slip effect, proved by both experiment [56][57][58] and molecular dynamic(MD) simulations 59 . Javadpour and Afsharpoor et al. [60][61][62] uncovered the influence of the slip effect on the infiltration of micronanopores.
A phase portrait is a geometric representation of the trajectories of a dynamical system in the phase plane. Phase portrait analysis is a remarkable tool for analysing the dynamic process of a complex system, such as spontaneous imbibition. The phase portrait method can evaluate the dynamic process of the system without solving the complex governing ordinary differential equations (ODEs). For mechanical systems, the phase space usually consists of all possible values of position and momentum variables. The concept of phase space was developed in the late 19th century by Ludwig Boltzmann, Henri Poincaré and Willard Gibbs 63 and is commonly used to analyse the possible state of a dynamical system.

Spontaneous Imbibition Model
Spontaneous imbibition occurs in capillaries as wetting fluid imbibed in capillaries. This process is driven by capillary force automatically. Liquid slippage in nanoscale capillaries should also be considered. Commonly, hydrophobic surfaces tend to cause water slip 52,53 . But hydrophobic nanopores will not lead imbibition, which is not the focus of our work. We just focused on the imbibition process in this work. In weak hydrophilic nanopores, water slip also occurs. Wu et al. 54,55 established a universal model for water flow through nanopores, effective slip is a linear sum of true slip and apparent slip. True slip is always positive. For water in hydrophobic tube, apparent slip is positive, while for water in hydrophilic tube, apparent slip is negative. Similar with Wu's model, in weak hydrophilic nanopores, wettability will give slight negative effect on effective slip, but not change effective slip into negative. This is what we discussed in this model. In strong hydrophilic tube, the tube wall will induce and locked layers of water molecule, which is a multilayer sticking model with no-slip boundary 54 . Our model for slip boundary is not suitable.
Also, contact angle changes with imbibition speed which makes the imbibition process additionally complex. This section shows the establishment of a spontaneous imbibition model while considering liquid slippage and contact angle variation.
Spontaneous imbibition only exists in wetting tube, only wetting phase can induce imbibition. For inorganic pores, imbibing wetting phase is water, while for organic pores, imbibing wetting phase is oil. To simplify the investigating object and focus on the relationship between the dynamic angle and slip length, we simplify the model as wetting phase liquid imbibed in weak wetting single tube with fractal tortuosity. Certain simplifications are made. a) The cross section of the tube is circular, b) the liquid is in the wetting phase, and gas is nonwetting, c) the liquid is Newtonian with laminar flow(Inertial forces are ignored), d) the driving force of spontaneous imbibition is capillary force and e) slip occurs at the tube wall, f) Tube wall is weak wetting for liquid. Figure 1 illustrates a schematic diagram of spontaneous imbibition in fractal tortuous capillaries. In the nanoscale system, slip length L s is comparable to capillary diameter. Hence, the slip effect should be considered. Velocity distribution on the cross section also follows a parabolic curve, assuming that a slip length L s exists. Liquid slippage can greatly increase liquid flux. On the basis of the Hagen-Poiseuille equation, fluid flux considering liquid slippage is as follows: where Q is fluid flux in the tube, Δp is the pressure difference of the fluid, μ is dynamic viscosity, L f is the length of tube based on liquid pathline, R is the equivalent radius of the tube and L s is the slip length. Oftentimes, hydraulic diameter is introduced to non-circular tubes as follows: c where λ is the hydraulic diameter, A c is the cross-sectional area and C is the wetted perimeter. For circular tubes, λ denotes the diameter or λ = 2R.
Δp represents the driving pressure. For upward flow, meniscus height increases during the imbibition process.
c 0 where L 0 is the distance between meniscus and tube intake. For straight and vertical tubes, L 0 is the height of meniscus, ρ is the liquid density, g is the gravitational acceleration, Ψ is the angle between the horizon and the tube and p c is the main driving force. According to the Young-Laplace equation, p c is as follows: c where σ is the interfacial tension between the gas and the liquid, θ is the contact angle and r is the tube radius. In non-circular tubes, p c is expressed as follows: During the imbibition process, the contact angle is dynamic. We select the liner relationship here. According to Hilpert 42 , contact angle θ in the Young-Laplace equation should be dynamic and expressed as follows: eq f where θ eq is the equilibrium contact angle and α ≥  0 is a non-dimensional parameter. Ca is a capillary number. This model, as Hilpert 42 assumed, is a special case of the power law and series models. Hence, p c for dynamic angle can be expressed as follows: c h eq f Liquid flux Q considering slippage, dynamic contact angle and gravity can be expressed as follows: Imbibition velocity is determined as follows:   Capillaries are not always straight. Thus, tortuosity should be considered. Yu and Cheng 20 introduced fractal dimension to express tortuosity and expressed the relationship between L f and L 0 as follows: where L f is the length of fluid pathline or real length of the capillary and L 0 is the linear distance between the meniscus and tube intake. Derived by Cai et al. 22 , the relationship between velocity v f and linear velocity v C can be expressed as follows: where λ is the capillary diameter, D T is the fractal dimension of a tortuous capillary and D T ranges from 1 to 3. For a straight capillary D T = 1, a higher D T represents higher tortuosity and longer pathline capillary. By combining Eqs (10) and (12), we can derive the ODE of imbibition as follows: (13) is a nonlinear implicit equation with independent variable t and dependent variables L 0 and =  L dL dt / 0 0 . The differential expression is as follows: At t = 0, the imbibition length is 0. Hence, the initial condition is as follows: By introducing dimensionless groups, the implicit ODE can be simplified as a dimensionless form to make the generalised equation applicable in different scales, which is helpful in performing sensitivity analyses. On the basis of the dimensionless method by Hilpert 42 , the three independent dimensionless groups are as follows: where L D is the dimensionless imbibition length, t D is the dimensionless time and L sD is the dimensionless slip length. On the basis of this dimensionless method, dimensionless velocity can be expressed as follows: We can non-dimensionalise the imbibition speed as follows: (18) remains an implicit ODE. By rearranging the first-order term, Eq. (18) can be expressed as follows: or simplified as follows:  (20) is an implicit first-order nonlinear ODE that represents the imbibition distance with respect to imbibition time considering liquid slip and dynamic contact angle. At initial condition t D = 0, the initial interface position is at the pressure intake. Hence, dimensionless imbibition length is expressed as follows: Analytical or numerical solutions are easy to obtain. In Eq. (18), the concept that the capillary force item should be always larger than gravity item should be stressed. When the capillary force is balanced against gravity, the meniscus reaches the equilibrium height, which is expressed as follows: Dend eq 2 Solutions Initial stage of imbibition or horizontal flow. In the case of horizontal flow Ψ = 0 or at the first stage of the imbibition process, capillary force in the nanoscale tube is much larger than that in the gravity item. The gravity item can thus be neglected. Hence, after ignoring the linear term ρ λ σ Ψ L g D sin 4 2 , the item can be simplified as follows: The equation is easily solved by separating the variables and subsequently implementing an integration method to derive an analytical solution. Taking the initial condition | = , the analytical solution is as follows: and solved by the separation of variables method, the analytical solution is as follows:   For other occasions, D T is not an integer. To the best of our knowledge, deriving the analytical solution is difficult. As the standard differential form has been derived as that in Eq. (19), the implicit ODE can be solved using a numerical method. We use the variable-step fourth-fifth-order Runge-Kutta algorithm (ODE45). This method is an adaptive-step numerical analysis process for ODEs and effective in obtaining the numerical solutions.

Summary of solutions.
Under different occurrences, we derived analytical and numerical solutions for spontaneous imbibition in fractal tortuous capillaries while considering dynamic contact angle and slip effect. Together with cases of gravity-ignored and -included processes and three categories of D T , we summarise the solutions. Table 1 illustrates the analytical and numerical solutions that we derived. Expressions of imbibition speed and length are included.

Phase Portrait Analysis
Phase portrait analysis is a significant tool for analysing the dynamic state of such a system as spontaneous imbibition. The main advantage is that the method can evaluate the system without solving the complex governing ODEs, which control complex dynamic systems. Obtaining explicit solutions that are expressible in finite terms is sometimes not easy for ODEs. A qualitative study has been conducted on differential equations in which the important characteristics of the solutions to differential equations were deduced without actually solving them 64 . In the following section, a geometrical method, the phase plane, is introduced to analyse solutions to the spontaneous imbibition process in fractal tortuous capillaries considering dynamic contact angle and slip effect.
A phase portrait is a geometric representation of the trajectories of a dynamical system in the phase plane, which consists of the system trajectories (arrows), stable steady states (dots) and unstable steady states (circles) in state space. The phase portrait method is useful in studying dynamical systems and forms a series of typical  Gravity included (analytical solution with implicit form) (analytical solution with implicit form) For spontaneous imbibition, we develop a series of phase portraits to analyse the stability of solutions to the imbibition model and influence of variables on the imbibition process. We first draw the phase line of spontaneous imbibition, as shown in Fig. 2(a), and then compare. We numerically solve Eq. (19) to derive the dimensionless imbibition length L D versus dimensionless imbibition time t D using the ODE45 method, as shown in Fig. 2(b). The negative value of L D does not have physical meanings as imbibition distance L D is always a positive value (L D > 0). Hence, negative values are ignored. Phase line is a 1D phase system that represents an interval of the domain of the derivative in the system. Values on the line denote imbibition distance, whereas arrows mean direction of the meniscus movement. As arrow density increases, velocity of meniscus decreases. Figuratively speaking, phase line Fig. 2(a) is a projection of Fig. 2(b). At the end of the imbibition speed, capillary force is equal to the gravity of the liquid. The dynamical system becomes stable, and this stage is also called the equilibrium state. , the gravity of the liquid is larger than the capillary force and the meniscus moves downwards (supposing the contact angle equations are the same). Finally, the dynamic system becomes stable at equilibrium height. The solution is drawn as a dashed line in Fig. 2(b) because solutions to L D that are larger than the equilibrium height are mainly not discussed in this model. If we set gravity as zero in this system, then the system becomes an unlimited one and L D will increase without limitation (if the tube length is infinite). Hence, the system will not reach a stable state. The solution line for an unlimited system is drawn as a blue line in Fig. 2(b).
A 2D phase figure composed by a pair of values  L L ( , ) D D is called a phase plane, which gives the imbibition velocity  L D at meniscus position L D . The solutions to the differential equation are a family of functions. Figuratively speaking, these solutions can be drawn in phase plane filled with a vector field. Vectors represent the direction of meniscus movement. Figure 3 shows a typical 2D phase plane and vector field for the solutions to Eq. (19). Given the initial condition, the solution line in phase plane is drawn as a red line. A certain part of the coordinate L D is omitted. The left vertical asymptote is L D = −D/C, and the right intersection with coordinate L D is L D = A/B. The origin of coordinates Q1 is a saddle point because vectors around point Q1 are opposite. At point Q2, where L D is larger than A/B, the velocity is a negative number. Hence, meniscus flow direction will flow back to Q2 and the dynamic system will be stable at Q2, which is a stable node. From the physical meaning of this model, the dashed line has no physical meaning, although it is a solution line. Area 0 < L D < A/B is our focus.

Model Verification
We use Hilpert's public data 41 of the spontaneous imbibition model to illustrate and verify our solutions for the imbibition process. Hilpert developed analytical solutions to spontaneous imbibition in capillaries without surface slip. By setting the same parameters, the solutions and published data are drawn under the same coordinate. Figure 4 shows the validation for dimensionless imbibition length and speed and dynamic contact angle.
Several physical properties of capillaries show their influence on imbibition height. Hilpert's model is a straight capillary without surface slip. By setting D T = 1 and L sD = 0, Eq. (19) can be simplified as Hilpert's model. Moreover, similar to Hilpert's model, α = .
 18 5, static contact angle θ eq = 70 0 . If contact angle remains stable in the

Results and Discussions
Slip length. In nanopores, the slip effect has a significant influence on the spontaneous imbibition process. Figure 5(a) shows the phase plane analysis under slip effect. When gravity is considered, the phase line extends to the X axis. By changing the dimensionless slip length of the system, a series of phase paths exists. All of the phase paths approach the X axis ( =  L 0 D ), which indicates that the dynamic system will be stable with any slip length. Furthermore, is an asymptote. Imbibition speed increases with slip length. Figure 5(b) shows the dimensionless imbibition length versus imbibition time for different dimensionless slip lengths (L sD ). As shown in this figure, the higher the L sD value, the longer the imbibition at the same imbibition time. Figure 5(b) shows that at t D = 2.5 × 10 7 , imbibition length L D is only 1700 for no-slip cases. The system reaches stable state at t D = 2 × 10 8 . For slip-influenced cases, as dimensionless slip length increases, average imbibition speed and dimensionless height increase. At t D = 2.5 × 10 7 , when L sD reaches 0.4, imbibition length L D becomes 2500, whereas L D only reaches 1700 for no-slip cases. Imbibition length L D increases by 47% because of slip length.
The simulated data indicate that the slip effect has an apparent advantage on the imbibition process. Most throats in shale are in microscale and nanoscale. When analysing the spontaneous imbibition process, overlooking slippage will considerably underestimate the imbibition flux, which is special compared with sandstones. Pore space in sandstones is large, with most diameters of throats in millimetre-scale and microscale, whereas slip length mainly exists in nanoscale 62 . Under this circumstance, the value of L sD is merely in a 0.001 scale, which is close to 0. Traditional imbibition models are based on imbibition analyses of sandstone, in which ignoring slippage is acceptable. Meanwhile, in imbibition analyses for shale, throats are mainly in nanoscale, which are in the same scale as slip length. Hence, dimensionless slip length L sD is approximately 0.1-10. Figure 5(b) shows that slip effect in this situation noticeably and significantly promotes the imbibition process. Hence, slippage cannot be neglected in micronanopores.
Dynamic contact angle. Apart from slippage, the dynamic contact angle has a significant influence on the imbibition process. Figure 6(a) shows the phase plane analysis of the influence of dynamic contact angle on spontaneous imbibition. When gravity is considered, the dynamic system can also reach equilibrium state because all of the phase paths point towards the equilibrium point. The dynamic contact angle is influenced by imbibition speed. Hence, capillary forces have a nonlinear relation to imbibition length. If α =  0, then the dynamic contact angle is equal to the static angle. As α increases, the effect of dynamic contact angle on imbibition process becomes larger. Figure 6(a) shows that dynamic contact angle mainly influences the initial stage of imbibition. The influence diminishes as imbibition length increases. Meanwhile, Fig. 6(b) shows the dimensionless imbibition length versus time for different dynamic coefficients (α). α significantly influences the early stage of imbibition. If t D > 40, then three curves are nearly parallel to one another, which indicates that their imbibition speeds are similar.
In the imbibition process, the rough flowing wall has resistance to fluid, which causes a 'delay' effect on the contact angle. The dynamic contact angle becomes bigger than the static one, whereas the capillary force in the dynamic process becomes smaller than that under static situation. High imbibition speed will cause a heavy 'delay' effect and small capillary force. At the initial stage, the imbibition speed or . L D is fast. Hence, the effect of the dynamic angle is apparent. At the end of the imbibition process, the effect of the dynamic angle is negligible. Gravity. In the previous sections, gravity is considered while the dynamic system is stable at equilibrium height. In a no-gravity environment, the dynamic system is not stable at the end. Discussing the stability of the system is easy using phase portrait analysis. Figure 7(a) shows the phase plane analysis with gravity. The phase paths that form a family of parallel straight lines are g = 0. The phase paths are apparently parallel straight lines that do not show any sign of convergence. Meanwhile, when gravity is considered, all phase lines converge at the stabilised line. Figure 7(b) shows the imbibition curve with and without gravity. The curve is solved using the ODE45 computation. As imbibition time increases, the gravity-considered imbibition height becomes stable and stops increasing at equilibrium height. Under the no-gravity condition, the increase in imbibition height does not stop. This result is the same as that in phase plane analysis.
Suggesting slip length is constant. The imbibition process varies with different capillary diameters. A large capillary gains a small L sD . The influence of gravity is noticeable in the imbibition process. At the beginning of the process, imbibition speed is much smaller in a large capillary than that in a small capillary regardless of whether gravity is considered or not. As L D increases, gravity can decrease the imbibition speed. In conclusion, gravity is one of the key factors that control the stability of this spontaneous imbibition system. The coupling relationship among gravity, dynamic contact angle and slip length is complex. To explicitly express this coupling relationship, 81 independent simulation cases that are based on the solutions to the imbibition process are performed in a 3D parameter space. Figure 8(a) and (b) show the parameter planes. Figure 8(a) depicts the parameter space at t D = 20, which indicates the initial stage of the imbibition process, and Fig. 8(b) shows the parameter space at t D = 2 × 10 8 , which indicates the last stage.
The figures show that at the beginning of imbibition, the two surfaces almost overlap and the influence of gravity is not large. α acts as a resistance, and L sD acts as a promotion. Meanwhile, a small α and large L sD characterise a positive imbibition process. The red line denotes a no-slip condition. Imbibition speed at the slip boundary is several times bigger than that at the no-slip boundary. At the end of the imbibition process, α and L sD have nearly no effect on imbibition height because the system has reached equilibrium height without the influence of these  parameters. This condition verifies the analytical solution previously discussed. Under the gravity-considered condition, α has no influence on imbibition height and L sD also acts as a promotion. Table 1, the spontaneous imbibition function can be derived from an analytical solution under certain special situations, such as D T = 1, 2, g = 0. Without these special situations, the imbibition function can be solved using numerical models, such as the ODE45 method, which is commonly used to solve complex ODEs. In Fig. 9(a), we compare the analytical and numerical solutions of the ODE for the imbibition process and obtain excellent agreement.

Typical situations and time complexity. As shown in
The fundamental process of the solution is iteration because ODE45 is a discretisation method. Although numerical solutions agree with the analytical solutions, running numerous iterations consumes considerable calculation time. Figure 9(b) shows the time consumed in one capillary for both solutions. Numerical solutions consume more time than analytical solutions. For certain complex ODEs, numerical methods do not easily converge during the iteration process and numerical solutions cannot be derived. For these ODEs, the phase portrait analysis method is a brilliant tool for system analysis because it does not rely on analytical solutions to complex ODEs. A large amount of calculation time can be avoided. Hence, phase portrait analysis can significantly shorten the calculation time and provide a promising way for complex dynamic system analysis.

Conclusions
In this work, the spontaneous imbibition model in fractal tortuous capillaries is successfully built, considering dynamic contact angle and slip effect during the imbibition process. On the basis of this model, ODEs are built describing dimensionless imbibition height versus dimensionless imbibition time. We introduce phase portrait analysis to analyse the imbibition process. Furthermore, when D T = 1, 2 and g = 0, analytical solutions for the imbibition equation are derived. Under other conditions, numerical solutions are derived. Additionally, in the imbibition model, we explore influential factors for spontaneous imbibition. The accuracy and reliability of phase  portrait analysis and the analytical solutions are examined. The phase portrait analysis method is a promising way for analysing complex spontaneous imbibition. Finally, a 3D parameter space based on the imbibition model is performed to express the relationship among slip length, dynamic angle and gravity.
Findings show that phase portrait analysis is a significant tool for analysing the spontaneous imbibition process. The most prominent advantage is that the phase portrait method can evaluate the process without solving the complex governing ODEs. Phase portrait analysis and the solutions to the ODEs reveal that dynamic contact angle and slip effect play an important role in fluid imbibition in fractal tortuous capillaries. The calculations show that considering slippage in imbibition simulation significantly increases imbibition flux in nanotubes, whereas considering dynamic contact angle corrects the imbibition process at the initial stage. Finally, gravity is a key factor that controls the stabilisation of the imbibition process.