Characteristic Sensitivity of Turbulent Flow within a Porous Medium under Initial Conditions

Flows within porous media play important roles in many scientific and industrial systems. However, the case wherein such flows become turbulent has not been completely understood, particularly, from a mathematical viewpoint. In this study, the $k-\varepsilon$ model (the variable $k$ denotes the turbulent kinetic energy per unit mass and $\varepsilon$ the dissipation rate for the turbulent kinetic energy), which has been widely used for the usual turbulent flow in a clear fluid, was applied to a turbulent flow through porous media. If a homogeneous, averaged flow and homogeneous isotropic turbulence are assumed, the governing equations for the variables $k$ and $\varepsilon$ describe a straight line as a nullcline common to both the variables. The nullcline was shown to be a line attractor. A temporal evolution of the eddy viscosity was also obtained, and the initial and final values of the eddy viscosity were observed to be related to a power law, indicating universal sensitivity. This sensitivity originates from the common nullcline and is not observed for the usual turbulent flow. Finally, nonlinear mathematical and seismological implications were provided using based on the results obtained.


I. INTRODUCTION
Porous media, or those media with a solid skeleton and many pores inside, are substantially studied in several scientific and industrial fields, including oil reservoirs, rockfill dams, fluidized bed combustion, biomaterials, filtration membranes, and chemical wastes [1][2][3][4][5][6]. For these aforementioned cases, the pores are considered to be filled with fluid, i.e., water and/or oil. A particular topic of research interest in these fields seems to be directed at the fluid flow and transport phenomena [7][8][9][10]. The most commonly adopted principle for many studies on these topics is the Darcy law, which states that flow rate through a porous medium is proportional to the fluid pressure gradient. This law assumes a laminar flow.
Porous media have also been of interest in seismological studies [11][12][13][14][15][16][17][18] mainly because fault rocks characteristically contain pores, and fluid (water) is considered to fill these pores. Particularly, an increase (or decrease) in the fluid pressure induces decrease (or increase) in the effective normal stress that acts on the fault zone, thereby decreasing (or increasing) friction stress and increasing (or decreasing) slip velocity. Many processes, such as thermal pressurization and dilatancy, are the presumed mechanisms causing fluid pressure variation [11][12][13][14][15][16][17][18]. Additionally, the increase (decrease) in the fluid pressure in the fault zone induces the fluid outflow to (inflow from) the outside of the fault zone because of the emergence of the gradient of the fluid pressure profile along the direction normal to the fault zone. This outflow (inflow) also reduces (raises) the fluid pressure in the fault zone, thereby affecting the slip behavior [14][15][16][17].
In fluid mechanics, turbulent flow in a clear fluid (the fluid in the system consisting of only the fluid phase, not the fluid in porous media) has attracted significant * t-suzuki@phys.aoyama.ac.jp research interest. In actuality, two main problems exist for the treatment of the turbulent flow. The first is the complexity involved in exactly simulating a real flow field because the degree of freedom for the flow is significantly high. The second is the sensitive dependence of the system behavior on the initial conditions, as widely known for turbulent flows; i.e., the Navier−Stokes equation has sensitive dependence on the initial conditions [19]. The viewpoint of time average or ensemble average has been introduced to circumvent the aforementioned complexity (i.e., the first problem) and sensitivity (i.e., the second problem). Many frameworks using such averages have been constructed; these frameworks include the mixing length model [20,21], k − ε model (where k denotes the turbulent kinetic energy per unit mass and ε the dissipation rate of k into a larger wave-number space) [22][23][24][25][26][27][28], and the Reynolds stress transport model [29][30][31]. Among these models, the k − ε model has been widely used. This model requires the decomposition of the fluid velocity field into the time-averaged velocity and deviation from it. The variables k and ε are expressed in terms of the deviation and used to calculate the eddy viscosity ν T , which is a quantitative measure of how momentum is transported from the time-averaged velocity field to the deviation thereof.
Turbulent flow in porous media has also been widely studied from theoretical  and experimental [55,56] viewpoints. Of the models that treat turbulent flow, the k − ε model has been widely employed [33,34,[36][37][38][39][40][42][43][44][45][46][47][48][49][50][51][52][53]. To consider the temporal evolutions of k and ε in porous media, it should be noted that energy is supplied from the time-averaged velocity to its deviation via pore walls because the averaged flow creates disturbance upon hitting the pore walls. Therefore, finite k and ε values in the steady state with constant time-averaged flow can exist in porous media [38], which is not the case for clear fluids. Such a statement implies that there can exist infinitely many steady states; i.e., attractors are continuously distributed on the k − ε phase space. Obtaining solution orbits connecting the initial states and attractors is a useful step to understanding the temporal evolution of the turbulent flow in porous media. Particularly, the sensitive dependence of the solution orbits on the initial conditions is expected by analogy with a previous nonlinear mathematical treatment [18], which revealed that a continuous attractor can produce sensitivity to the initial conditions.
Previous studies assumed variations in the initial conditions for k and ε in porous media. For example, several researchers assumed the value k 0 /ε 0 (where k 0 and ε 0 denote the initial values of k and ε, respectively) to be constant [23,25]. Another study [39] employed the "constant eddy viscosity" condition that was initially introduced for clear fluids [27]. Lin and Liu [27] introduced a small amount of k 0 in the initial condition as a seed for generating turbulence in clear fluids, whereas ε 0 was determined by the condition that the eddy viscosity was constant. We should systematically treat these conditions.
Continuous attractors have drawn the interest of many researchers, especially those from the field of neural networks [57][58][59][60][61]. For example, line attractors have been explored for studying oculomotor control. However, the mathematical understanding of such systems, particularly concerning the effect of the disturbance in the initial conditions on the steady state, has not been performed thus far. To understand this effect, the analogy of a turbulent flow in porous media works well. This paper is organized as follows. Our model and its governing equation system are presented in Sect. II. A homogeneous, time-averaged flow and a homogeneous, isotropic disturbance are assumed, and the analytical treatment for the state is discussed in Sect. III. In this section, the universal sensitivity of the final state to the initial state is emphasized for the turbulent flow in porous media. Moreover, the sensitivity of the final eddy viscosity on the initial eddy viscosity, and that of the difference in the final eddy viscosity on the difference in the initial eddy viscosity are clarified. The findings in this study, as well as some nonlinear mathematical and seismological implications, are summarized in Sect. IV.

II. MATHEMATICAL FORMULATION
We now present a mathematical framework that governs the behavior of turbulent fluid flow in a porous medium. Let us consider an infinite porous medium, whose pores are filled with fluid, i.e., water. The solid phase of the medium is assumed to be rigid such that the porosity φ remains constant. The fluid phase is assumed to be single-phase Newtonian fluid and incompressible such that the density of the fluid phase ρ remains constant. Moreover, consider the pore Reynolds number Re p , expressed as Re p ≡ dv c /ν, where d, v c , and ν denote the characteristic pore scale, characteristic flow velocity scale, and viscosity of the fluid phase, re-spectively. Reportedly [38,62], a flow can be considered turbulent when where Re c denotes the critical pore Reynolds number, above which the flow becomes turbulent. Generally, Re c is considered to be in the range of 300 − 1000 [34,47,62].
To deal with such turbulent flow within porous media, we utilize the k − ε model constructed by [33,38]. First, let us suppose that the length scale is smaller than the characteristic pore scale, i.e., from a microscopic viewpoint. The framework treating the "usual turbulent flow," or the "turbulent flow in a system consisting of only the fluid phase (clear fluid)" can be used for this length scale. We write k (turbulent kinetic energy) and ε (dissipation rate of k into a large wave-number space) for the usual turbulent flow as k usual and ε usual , respectively, and the mathematical framework of the k − ε model in porous media is constructed from k usual and ε usual in the following discussion.
We then decompose the flow velocity profile into timeaveraged velocity and deviations from it, as mentioned in Sect. I. The time average of variable A, A, is defined as where ∆t denotes the integration time interval. Performing this integration eliminates the high-frequency components of the variable. Using Eq. (2), the velocity field u is expressed as where u ′ is the deviation of u from u. Notably, u ′ = 0 by definition. For the usual turbulent flow, k usual and ε usual can be written in terms of u ′ as [23,29] and where i, j = 1, 2, 3 and the repeated indices are summed. Equation (4) is based on the following assumption: where ν T denotes the eddy viscosity and δ ij the Kronecker delta. The term u ′ u ′ is known as the Reynolds stress, and Eq. (6) shows that the momentum is transported from the time-averaged velocity u to its deviation u ′ . Equation (6) also states that the Reynolds stress is directly proportional to the local time-averaged velocity gradients and that the proportionality factor, i.e., the eddy viscosity ν T , is a scalar quantity. Notably, ν T quantitatively characterizes how the flow is disturbed and is regarded herein as the intensity of disturbance.
Nonetheless, we must consider a solid skeleton while discussing k and ε in porous media. Thus, we must consider a coarse graining viewpoint associated with the representative elementary volume (REV), as shown in Fig. 1. The volume of REV, ∆V , is several times larger than the characteristic pore scale and significantly smaller than the system size; thus, it becomes necessary to take the spatial averages of the variables within the REV. Accordingly, the spatial average of variable A, A , is expressed as where ∆V f denotes the volume of the fluid phase within ∆V . We, therefore, conclude that the temporal evolutions of k usual and ε usual must be obtained to understand the behavior of the turbulent flow within porous media.
The authors of [38] derived the governing equations for k usual and ε usual . Their framework, which was used by other studies [43,46], is also employed in our study as follows: where K denotes the permeability of the porous medium, , and C ε2 are model parameters. The velocity field u D and permeability are related as which is the well-known Darcy law, with p denoting the fluid pressure. Hereinafter, we will refer to u D as the Darcy flow and the velocity field u ′ as turbulence. Finally, we must define the eddy viscosity ν T to close the macroscopic governing equation system [38] as where C µ is an empirical constant. The first terms in the right-hand sides of Eqs. (8) and (9) describe the diffusion effects of k usual and ε usual , respectively. The second terms, including the Reynolds stress, are production terms, which express the energy conveyed from the Darcy flow into the turbulence. The terms −ρφ ε usual and −C ε2 ρφ ε usual 2 / k usual denote the energy dissipation due to the fluid viscosity. Such terms also appear in the governing equations for the usual turbulent flow with the limit φ → 1. Additionally, other production terms in porous media are introduced, such as the terms including K in Eqs. (8) and (9). Notably, 1/ √ K is a measure of the flow resistance [see Eq. (10)], and the higher kinetic energy of u D will be converted to k usual with higher 1/ √ K via the pore walls. Noteworthily, the usual turbulent flow vanishes at the steady state for homogeneous u D because the term with 1/ √ K does not exist. For porous media, the turbulence can continue to remain in such a case because the energy is always supplied from the Darcy flow into the turbulence via the pore walls, which plays an important role in this study.
For the usual turbulent flow, the parameter values are widely assumed as C µ = 0.09, σ k = 1.0, C ε1 = 1.44, C ε2 = 1.9 − 1.92, and σ ε = 1.2 − 1.3 [22,25,28] and are considered universal, i.e., they do not depend on the flow. We will use these usual parameter values (C ε2 and σ ε are assumed to be 1.9 and 1.2, respectively) assuming that the flow is close to the usual turbulent flow as the first approximation [36,47]. Additionally, we assume c k = 0.28, which was numerically derived by solving the Navier−Stokes equation for the flow through arrays of circular rods and by imposing the condition that the obtained results can be reproduced by using the macroscopic equations for k usual and ε usual [38,44,51].
However, obtaining the value of c k is somewhat difficult, and thus some researchers insist other values, e.g., c k = 0.16 − 0.34 [63], and c k = 0.22 − 0.34 [44]. These differences are based on the porosity and configurations of pores. Nonetheless, the following discussion does not qualitatively change while choosing values in the range suggested for c k by previous studies. This is because the solid phase is assumed to be rigid and the porosity and the pore configurations remain unchanged, thereby indicating that parameter c k does not change temporally.
Finally, we consider that k usual and ε usual correspond to k and ε in porous media, respectively. Therefore, we will express k usual and ε usual as k and ε, respectively, in the following discussion. Equations (8) and (9) are considered to be the governing equations for k and ε.

A. Nullclines and Solution Orbit
We assume a homogeneous isotropic porous medium, and the constant Darcy flow, which is defined as a flow whose velocity is constant everywhere. This assumption has been considered for both the usual turbulent flow [23,24,26] and turbulent flow in porous media [33,38,63]. For example, Zhao [26] insisted that such a situation is realized in the region outside the boundary layer over a flat plate, as shown in Fig. 2. Here, we choose the x axis along the Darcy flow direction. We can consider that a constant fluid pressure gradient is maintained along the x axis by providing fluid from an infinitely far point. It is also important that u D can be written as u D = u Dx , where u D denotes a positive constant andx a unit vector directed along the positive x axis. The turbulence is assumed to be homogeneous and isotropic. We also assume ρ and φ to be constant, as mentioned in Sect. II.
Under this condition, we can neglect the spatial differentiations in governing Eqs. (8) and (9); therefore, we can express the governing equations for k and ε as respectively. We should emphasize that exactly the same straight line given by the equation, is a nullcline common to both the variables on the k − ε phase space. Hence, the solution orbits are neither horizontal nor vertical where the solution orbits and this nullcline intersect [18]. Moreover, this line is a line attractor or repeller, as observed in other systems including [18]. Furthermore, the straight line ε = 0 (k axis) is a nullcline for ε.
We now obtain the analytical forms of the solution orbits. Equations (12) and (13) yield We assume constant k 0 = 0 and ε 0 = 0, where k 0 and ε 0 are the initial values of k and ε, respectively. Accordingly, Eq. (15) easily leads to the solution as follows: We now show that the common nullcline ε = c k u D k/ √ K is a line attractor. Figure 3 shows the direction of the temporal evolutions of k and ε. First, we define Region I as the region 0 < ε < c k u D k/ √ K and Region II as the region ε > c k u D k/ √ K > 0, in the k − ε phase space. We also define (k f , ε f ) as the point where the solution orbit (16) intersects the nullcline. With these definitions in hand, we obtain the relations k 0 < k f and ε 0 < ε f if (k 0 , ε 0 ) is in Region I because ε ∝ k on the nullcline, ε ∝ k Cε2 on the solution orbit, and C ε2 > 1. We can also conclude that k 0 > k f and ε 0 > ε f if (k 0 , ε 0 ) is in Region II. Second, we emphasize that k and ε increase (decrease) with increasing time in Region I (II), as the right-hand sides of Eqs. (12) and (13) are positive (negative). Therefore, if (k 0 , ε 0 ) is in Region I (II), the solution moves to the upper right (lower left) and gets absorbed into the point (k f , ε f ) on the nullcline with the limit t → ∞ (see Fig. 3). Thus, we can conclude that the nullcline is a line attractor, not a repeller. The steady stable solution is given by (k, ε) = (k f , ε f ), and this state will be referred to as the final state in the following discussions. Noteworthily, k and ε vanish with the limit t → ∞ for the usual turbulent flow, even though solution orbit (16) is valid for such flows [26]. Actually, the usual turbulent flow is described by the limit K → ∞. The nullcline is the k axis in such a case, with Region I vanishing. Therefore, all the solutions are absorbed into the origin. The terms including the finite K enable the homogeneous isotropic turbulence to survive as t → ∞ under the constant Darcy flow. Therefore, the relationship ε f = c k u D k f / √ K must be satisfied for the constant Darcy flow and fully developed homogeneous isotropic turbulence. Using this relationship and Eq. (16), k f and ε f are given as and respectively. These values are important in investigating the temporal evolution of the eddy viscosity, as discussed in Sect. III B.

B. Temporal Evolution of Eddy Viscosity
We now investigate the temporal evolution of the eddy viscosity. We note that the initial value of ν T , ν T 0 , is clearly given by and the value of ν T at the final state, ν T f , can be expressed as .
(20) This value is derived from Eqs. (17) and (18). The ratio between the final and initial viscosities is given by where α ≡ ε 0 /k 0 . Here, we emphasize −(C ε2 − 2)/(C ε2 − 1) > 0, as C ε2 = 1.9. Let us assume that (k 0 , ε 0 ) is in Region I. In this case, c k u D / √ K > α, and hence ν T f > ν T 0 . Alternatively, if we assume that (k 0 , ε 0 ) is in Region II, then ν T f < ν T 0 . Briefly, if (k 0 , ε 0 ) is in Region I, the final disturbance is stronger than the initial one, and if (k 0 , ε 0 ) is in Region II, the final disturbance is weaker than the initial one. Notably, ν T is a quantitative measure of the disturbance intensity, as mentioned in Sect. III A.
From a physical viewpoint, the energy supply from the Darcy flow to the turbulence is dominant in Region I, and energy dissipation is dominant in Region II. To elucidate on this, let us consider the straight line k = const. On the line, ε is greater in Region II than in Region I. Moreover, as ε denotes the energy dissipation, its effect dominates and consequently the disturbance is reduced in Region II.
The investigation above can be easily understood on the basis of k − ε phase space. We now investigate the temporal evolution of the eddy viscosity. We begin by obtaining the analytical form of the curve on which the eddy viscosity takes the same value (iso-eddy-viscosity curve) on the k − ε phase space. From Eq. (11), we have for the curve. This equation insists that the iso-eddyviscosity curve is a parabola whose vertex lies at the origin in the k − ε phase space. Comparing Eqs. (16) and (22), we can show that if (k 0 , ε 0 ) is in Region I, then the eddy viscosity increases with time. This increase in the eddy viscosity occurs because the gradient of the solution orbit is smaller than that of curve (22) at their crossing point. Notably, C ε2 = 1.9 < 2 (see Fig. 4). With increase in time, the solution migrates to the isoeddy-viscosity curve with a higher eddy viscosity than before. For the same reason, if (k 0 , ε 0 ) is in Region II, ν T decreases with time (see Fig. 4).

C. Relationship between k0 and ε0 and Universal Power Law
We allowed any values for k 0 and ε 0 in the above discussion because, particularly, determining ε 0 is an arduous task. However, several previous studies assumed that the values of k 0 and ε 0 are ruled by some simple relations. For example, some researchers assumed that the ratio ν T /ν is constant for the usual turbulent flow [27,64] and for the turbulent flow in porous media [39], leading to the relation ε 0 ∝ k 2 0 . Another example ε 0 ∝ k 0 was used for the usual turbulent flow [23] and for the turbulent flow in porous media [31]. The relation ε 0 ∝ k 3/2 0 was adopted for the usual turbulent flow [65]. Although several ideas were considered, we, at least, require "the larger ε 0 with larger k 0 ." Accordingly, we assume ε 0 ∝ k n 0 , where n denotes a positive number that satisfies the condition 1 ≤ n ≤ 2.
We first analytically show the dependence of ν T 0 and ν T f on k 0 . We define β as the coefficient of proportionality for ε 0 and k n 0 , i.e., ε 0 = βk n 0 . From this definition and Eqs. (19) and (20), we obtain the relations and Furthermore, if n = 1, we have the fixed point (k fix 0,n , ε fix 0,n ) on the k − ε phase space. This point does not move with increase in time. The point (k fix 0,n , ε fix 0,n ) is the point of intersection of the nullcline and the curve that describes the relationship between k 0 and ε 0 . From this definition, k fix 0,n is given by the solution of the equation which reads We can also obtain . (27) With these values, the fixed eddy viscosity, ν fix T,n , is given as .
(28) We now investigate the relationship between the normalized valuesν T 0 ≡ ν T 0 /ν fix T,n andν T f ≡ ν T f /ν fix T,n on the basis of Eqs. (23), (24), and (28) with n = 1 and 2. From these equations, we can obtain the relatioñ Using C ε2 = 1.9 and 1 < n < 2, we can prove that (C ε2 − n)/[(2 − n)(C ε2 − 1)] is positive when 1 < n < C ε2 . Therefore, for 1 < n < C ε2 ,ν T f monotonically increases with increasingν T 0 ; i.e., the stronger initial disturbance leads to the stronger final disturbance when 1 < n < C ε2 . On the other hand, for C ε2 < n < 2,ν T f monotonically decreases with increasingν T 0 ; i.e., the stronger initial disturbance leads to the weaker final disturbance when C ε2 < n < 2. If n = C ε2 ,ν T f is unity (ν T f = ν fix T,n ), which is independent ofν T 0 . With the k − ε phase space, we next show the sensitive dependence ofν T f on the initial conditions for the cases C ε2 < n < 2 and 1 < n < C ε2 . For the cases n = 1 and n = 2, we investigate such relations by employing Eqs. (23) and (24). We visualize the analytical result by using the k − ε phase space for the case C ε2 < n < 2. The solution orbits and the curve that describes the relationship ε 0 = βk n 0 are shown in Fig. 5. Notably, on the nullcline and red curve, ν T is higher on the points farther from the origin [see Eq. (22)]. Additionally, a point nearer to the origin on the red curve evolves to a point farther from the origin on the nullcline with the limit t → ∞ because the power of the solution orbit (C ε2 ) is smaller than n. We can, thus, conclude that a lowerν T 0 results in a higherν T f . The valueν T 0 = 1 is critical and fixed; i.e., ifν T 0 = 1, thenν T f = 1.
Moreover, the power law (29) implies the sensitivity of the system behavior within C ε2 < n < 2. Let us assume that (k 0 , ε 0 ) is in Region I andν T 0 ≪ 1. Here, note that ν T f is described by the negative power ofν T 0 , and that a significantly highν T f is induced by such a negligible value ofν T 0 . This sensitivity of the final state to the initial condition is peculiar to a turbulent flow in porous media. Although the usual turbulent flow is also sensitive to the initial condition, the sensitivity characterized by the power law (29) is specific to the porous media. This law is universal because the power (C ε2 −n)/[(2−n)(C ε2 −1)] does not depend on the media details including porosity and permeability.
We then consider the case 1 < n < C ε2 . We have visualized the behavior of the solution for this case in Fig.  6, where the solution orbits and the curve that describes the relationship ε 0 = βk n 0 are illustrated. On the nullcline and red curve, as previously mentioned, ν T is higher on the point farther from the origin. Additionally, a point nearer to the origin on the red curve evolves to a point nearer to the origin on the nullcline with the limit t → ∞ because the gradient of the solution orbit is larger than that of the red curve where both curves intersect (Notably, n < C ε2 ). In such a case, a higherν T 0 results in a higherν T f . Nonetheless,ν T 0 =ν T f = 1 remains a fixed value, as observed for the case C ε2 < n < 2.
Additionally, we can show that the universal sensitivity reappears. Because the power (C ε2 −n)/[(2−n)(C ε2 −1)] is positive and less than unity for 1 < n < C ε2 ,ν T f can be significantly higher thanν T 0 withν T 0 ≪ 1. Nonetheless, note that the sensitivities for the cases C ε2 < n < 2 and 1 < n < C ε2 are qualitatively different from each other. To clarify this, let us consider the limitν T 0 → 0. With this limit, we haveν T f → ∞ for the case C ε2 < n < 2 andν T f → 0 for the case 1 < n < C ε2 . We can insist that the sensitivity discontinuously changes for the cases C ε2 < n < 2 and 1 < n < C ε2 , and that the sensitivity is stronger in the former case than in the latter one. Finally, the common nullcline (14) is important for the universal sensitivity. The sensitivity due to the common nullcline has been shown to emerge in a dynamic earthquake slip process, as detailed in [18].
If n = 1, we cannot define the fixed point and cannot employ the power law (29). In this case, we have a relationship ν T f ∝ ν T 0 from Eqs. (23) and (24). The sensitivity does not emerge when n = 1.
If n = 2, we can confirm that ν T 0 has constant value C µ /β, and thus we cannot use the power law (29). Alternatively, we can employ Eq. (24), which represents the power law, as We can conclude that ν T f is sensitive to k 0 or ε 0 because (C ε2 − 2)/(C ε2 − 1) is negative. Several previous calculations [38,63] assumed a homogeneous Darcy flow and constant k 0 and ε 0 , and maintained their values at the boundary of a one-dimensional computational domain with increasing time. In these calculations, the relationship ν T 0 > ν T f was demonstrated wherein the distance from the boundary was so large that the assumption of the homogeneous isotropic turbulence held valid. However, the relationship ν T 0 < ν T f holds true in laboratory experiments such as the one in [56]. A finite ν T f is generated from the negligible ν T 0 value in the experiments. This variation observed in previous studies may be explained using the single framework constructed in this study.

D. General Treatment of Difference in the Initial Conditions and Its Application to Nonlinear Mathematics
If the condition ε 0 ∝ k n 0 is not assumed, we cannot derive relations (23) through (29). However, let us consider two points (k 0 , ε 0 ) and (k 0 +δk 0 , ε 0 +δε 0 ), where δk 0 and δε 0 denote positive or negative amounts that satisfy the conditions |δk 0 /k 0 | ≪ 1 and |δε 0 /ε 0 | ≪ 1, respectively, on the k − ε phase space. Actually, ν T f can be sensitive to δk 0 and δε 0 , and the sensitivity should be analytically clarified. Note that the sensitivity investigated above was that of ν T f on ν T 0 , while the sensitivity studied here is that of the difference in ν T f on the difference in ν T 0 .
To clarify the sensitivity, we consider the following two relationships: and Therefore, if we neglect the terms of orders δk 2 0 and δε 2 0 and higher, we obtain where g ≡ δε 0 /δk 0 . We refer to 1 + ε 0 (2 − C ε2 )/(k 0 g − 2ε 0 ) as γ, which is plotted in Fig. 7 as a function of g. If g is near 2ε 0 /k 0 , the slight disturbance in the initial condition is enlarged in ν T f , thereby indicating strong sensitivity. Note that γ is not a power, but the relative magnitude of the difference in ν T f to the difference in ν T 0 . Therefore, the sign of γ does not determine the strength of the sensitivity, and thus only the absolute value of γ is relevant.

IV. DISCUSSION AND CONCLUSIONS
For the turbulent flow through porous media, we applied the k − ε model and assumed the constant Darcy flow and the homogeneous isotropic turbulence. The eddy viscosities of the initial and final states are related via a power law, which is universal and independent of the media details including porosity and permeability. Particularly, the sign of the power depends on the initial condition; if ε 0 ∝ k n 0 is satisfied, it is negative (positive) for the case C ε2 < n < 2 (1 < n < C ε2 ). For both the cases, the final eddy viscosity is sensitive to its initial value and the sensitivity is stronger for the C ε2 < n < 1 than for 1 < n < C ε2 cases. Furthermore, even if n = 2, the final eddy viscosity sensitively depends on the initial values k 0 and ε 0 . The sensitivity mentioned here is not observed for the usual turbulent flow, but it originates from the common nullcline for k and ε, as consistent with a previous study [18]. The importance of the common nullcline should be emphasized in nonlinear mathematics. Moreover, the difference observed in the final eddy viscosity is also sensitive to that observed in the initial eddy viscosity.
Let us now consider a spatially inhomogeneous flow. In this case, the spatial differentiations in Eqs. (8) and (9) do not vanish. In particular, the coefficients of the diffusion terms include ν T , which enhances the diffusions of k and ε. Although this diffusion term is not considered in this study, we qualitatively employ the effect of disturbance enhancing the time-averaged velocity in this section as the first step toward extending the present results by considering seismological slip as an example. Quantitatively considering this effect will be an important future work.
As mentioned in Sect. I, fault rocks can be considered porous media. Two semi-infinite porous media in contact will slide when they experience shear stress, thereby mimicking the fault motion. Moreover, the contact surface can be regarded as a fault zone. Particularly, the contact surface cannot be a mirror-like plane and is rough. Therefore, here we use the term "fault zone" instead of the term "fault plane." Previous studies were able to model two effects, namely, thermal pressurization and dilatancy effects (both of which change the fluid pressure on the fault zone, p F ), into a single framework [13,[15][16][17][18]. Herein, we assume that p F is initially the same as the ambient fluid pressure. The thermal pressurization effect increases p F , whereas the dilatancy effect reduces p F . Moreover, p F and the slip velocity on the fault zone assume a linear relationship: an increase (decrease) in p F corresponds to an increase (decrease) in the slip velocity (for a detailed discussion on this, kindly see [16][17][18]). Moreover, the fluid flow also affects p F and induces slip velocity change, as mentioned in Sect. I.
Let us provide a rough estimation of Re p during a dynamic earthquake slip process. Importantly, ν = 3×10 −7 [m 2 /s] for heated (near 373.15 [K]) water, and we can roughly estimate d to be in the range of 0.1 − 1 [mm], as suggested by gouge sizes for natural fault rocks [66,67]. We should also evaluate the characteristic flow velocity v, which is arduous to calculate. Specifically, to evaluate v, we must introduce the fault-zone thickness T F and characteristic time scale (rise time) D and consider the value T F /D as a measure of v. For subduction thrust faults, the order of mm−cm seems reasonable for T F [68][69][70][71]. Noteworthily, D is on the order of 0.1 − 10 [s] for earthquakes with a moment magnitude scale in the range of 6 − 7 [72][73][74]. Therefore, we can roughly estimate v as T F /D ∼ 0.01 − 10 [cm/s]. The corresponding pore Reynolds number at these ranges is approximately 0.033 − 3.3 × 10 2 , with the latter satisfying condition (1). Accordingly, we can conclude that the turbulent flow plays a more important role for earthquakes with larger gouges, thicker fault zones, and shorter rise times.
The increase in ν T enhances the fluid inflow or outflow because effective viscosities, ν + ν T /σ k and ν + ν T /σ k , include ν T , as shown in Eqs. (8) and (9). For ordinary earthquakes, the fluid inflow is considered to occur [17], and if the inflow is enhanced, the slip velocity would increase further because the subsequent thermal pressurization would raise p F , indicating decrease in the friction stress. Therefore, the acceleration is considered positive. The acceleration of the initial phase of the slip may be explained by the feedback proposed herein.
We finally extend the result obtained in Sect. III D from the viewpoint of the nonlinear mathematics. We consider two quantities, x 1 and x 2 . We define a single scalar value V (x 1 , x 2 ) on the x 1 − x 2 phase space, and V is assumed to be proportional to x a 1 x b 2 , where a and b denote constant numbers. The initial values of x 1 and x 2 are denoted by x 10 and x 20 , respectively. The common nullcline describing the continuous attractor is assumed to exist on the x 1 − x 2 phase space. Actually, a continuous attractor is a straight line that crosses the origin in many cases [18,61], and this assumption is also employed here. Therefore, the relationship x 2f ∝ x 1f is satisfied. Additionally, in many cases, the solution orbits can be described by x 2 ∝ x N 1 , where N denotes a constant, and we also adopt this relation. With these assumptions, we now expect V 0 ≡ V (x 10 , x 20 ) ∝ x M0 10 x N0 20 and where g ≡ δx 20 /δx 10 . If we define γ as (N f x 10 g + M f x 20 )/(N 0 x 10 g + M 0 x 20 ) in this case, the value of γ determines the sensitivity of V on the initial condition. Continuous attractors have drawn the attention of researchers in terms of, particularly, neural networks [57][58][59][60][61]. For example, x i described the activity of neuron i in the model of [61]. When we consider two neurons, the straight line x 2 = x 1 is the common nullcline, and this case is known to generate N = −1 [61]. In their model, we can consider x 2 as an example of V . We then have V 0 = x 20 and V f = √ x 10 x 20 . Therefore, we have M 0 = 0, N 0 = 1, M f = 1/2, and N f = 1/2 in this case, leading to γ = (gx 10 + x 20 )/2x 10 g. Therefore, if g is negligibly small, a negligible but finite difference in the initial condition is enlarged in the final value. However, if g is near −x 20 /x 10 , the enlargement is not observed. The parameter γ can be a measure for the sensitivity of physical quantities on the initial conditions in a system with common nullclines.