Well-posed continuum equations for granular flow with compressibility and μ(I)-rheology

Continuum modelling of granular flow has been plagued with the issue of ill-posed dynamic equations for a long time. Equations for incompressible, two-dimensional flow based on the Coulomb friction law are ill-posed regardless of the deformation, whereas the rate-dependent μ(I)-rheology is ill-posed when the non-dimensional inertial number I is too high or too low. Here, incorporating ideas from critical-state soil mechanics, we derive conditions for well-posedness of partial differential equations that combine compressibility with I-dependent rheology. When the I-dependence comes from a specific friction coefficient μ(I), our results show that, with compressibility, the equations are well-posed for all deformation rates provided that μ(I) satisfies certain minimal, physically natural, inequalities.

and temporally dependent. These are governed by conservation laws plus constitutive relations. Conservation of mass gives the scalar equation (∂ t + u j ∂ j )φ + φ div u = 0, (2.1) and conservation of momentum gives the vector equation where ρ * is the constant intrinsic grain density and g is the acceleration due to gravity. Closure of these equations is achieved through three constitutive relations.

(a) The Coulomb constitutive model
For a Coulomb material, which is assumed to be incompressible, the first constitutive relation states that φ is a constant. This then reduces (2.1) to the Flow rule: div u = 0. (2.3) For the next constitutive relation, the stress tensor is decomposed into a pressure term (where p = −σ ii /2) plus a trace-free tensor τ called the deviatoric stress. The second relation is then the Yield condition: τ = μp, (2.5) where μ is a constant and for any tensor T the norm is defined by This yield condition expresses the idea that a granular material cannot deform unless the shear stress is sufficient to overcome friction. 1 The third constitutive relation requires that the eigenvectors of the deviatoric stress tensor and the deviatoric strain-rate tensor 2 D ij = 1 2 (∂ j u i + ∂ i u j ) − 1 2 (div u)δ ij (2.7) are aligned (see figure 1 for motivation), which may be written Alignment: In words, the above equation may be interpreted as asserting that in the space of trace-free symmetric 2 × 2 matrices, which is two-dimensional, D and τ are parallel. Thus, this matrix equation entails only one scalar relation. For reference below we record that It is customary [2,5] to process these equations by expressing the deviatoric stress τ in terms of p and the strain rate as follows: where we have invoked (2.5) and (2.8). We may substitute (2.10) into (2.2) to obtain (2.11) and the resulting equation, together with (2.3), gives three equations for pressure p and velocity u.
In form, at least, these equations resemble the incompressible Navier-Stokes equation. However, in two dimensions (as considered here) they are always ill-posed [2].

(b) Incompressible μ(I)-rheology
Work described by the Groupement De Recherche Milieux Divisés [4] has significantly improved the Coulomb model by including some rate dependence (in the sense of plasticity [29]) in the yield condition while making no changes in the incompressible flow rule (2.3) and the alignment condition (2.8). Specifically, a wide range of experiments are captured by replacing the constant μ in (2.5) by an increasing function μ(I) of the inertial number, (2.12) where d is the particle diameter. The expression μ(I) = μ 1 + μ 2 − μ 1 I 0 /I + 1 , (2.13) where μ 1 , μ 2 and I 0 are constants with μ 2 > μ 1 , is a frequently used form [30]. Below we shall assume that μ (I) > 0 and μ (I) < 0. (2.14) The modified yield condition changes (2.11) to read For the specific μ(I)-curve (2.13) this inequality covers a significant range of inertial numbers, specifically when the deformation rate is neither too small nor too large relative to the pressure. Outside of this range, the maximal-order linear stability analysis and numerical simulations show that perturbations grow exponentially with growth rates tending to infinity as their wavelength is reduced [9]. This behaviour is the hallmark of ill-posedness and leads to unphysical numerical solutions that strongly depend on the grid resolution used.

(c) Compressibility and I-dependent rheology
We refer to CSSM (see appendix A) for guidance in introducing compressibility into the rheology. Thus, we make no change in the alignment condition (2.8); we assume φ-dependence in the yield condition, and satisfy the inequalities and We may now state our main result, the well-posedness theorem for the system (

(d) Derivation of evolution equations
To place the equations in a larger continuum-mechanics context, we show that the CIDR equations of motion can be rewritten as a system of three evolution equations for the velocity u and the solids fraction φ. In form, these equations are analogous to the Navier-Stokes equations for a viscous, compressible fluid. We make no use of this form of the equations in our proof of wellposedness.
We want to eliminate stresses from the equations of motion. To this end, we propose to solve for the mean stress p using the flow rule (2.18), which we rewrite as Note that f (p, φ, I) depends on p both directly in its first argument and indirectly through I = 2d D / p/ρ * in its third argument. However, which by assumption (2.20b) is non-zero. Thus, we may apply the implicit function theorem to (2.21) to solve p = P(∇u, φ). 4 Given this, we may define and substitute into conservation of momentum to obtain an equation This equation, along with (2.1), gives a system of three evolution equations for the velocity u and the solids fraction φ. It is possible that previous formulations of compressible μ(I) equations [15][16][17] may be seen as CIDR equations with specific constitutive laws specified. In this paper, we choose to elucidate the well-posedness result with more generic choices of f and Y in order to remain impartial.

(e) Connection to μ(I)-rheology
Without making any attempt to be general, we illustrate one example of how μ(I)-rheology may be included in constitutive relations of the form (2.17), (2.18). Motivated by equation (A 3) in appendix A, we make the ansatz In these equations, it is worth emphasizing that p, φ, I are treated as independent variables, not to be confused with the dependence of I on p in the previous subsection. The function C(φ) is an increasing function of φ. As φ varies (with I fixed) the yield loci τ = Y(p, φ, I) derived from (2.24a) form a nested family of convex curves in stress space (figure 2b  (2.28) Solving this linear equation for α(I), with an integrating factor, we obtain  By (2.14), μ (I) < 0, so μ (sI) > μ (I) for 0 < s < 1. Thus, the last inequality using (2.14). Consequently, as desired.

(f) The incompressible limit
Based on an analogy with CSSM, let us suppose that C(φ) is a sensitive function of φ, say of the form where g is the constant of gravitational acceleration, b is a non-dimensional parameter and the factor ρ gd gives C the dimensions of pressure. The non-dimensional functionČ has an argument that is dependent on the minimum solids fraction φ min for sustained stress transmission between grains (random loose packing 5 ) and φ max , which is the maximum packing fraction that can be attained. Typically, φ = φ max − φ min is small. For definiteness we may takě In physical terms, the maximum solids fraction φ max represents the jamming threshold. We call the limit φ → 0 incompressible because, as may be seen from (2.32), the density of the material becomes essentially constant.
Next we use the yield condition to eliminate φ, reducing this number to four. Specifically, substituting (2.24a) into (2.17), we write the yield condition where the dependence on ∇u comes from the fact that I = 2d D / p/ρ * . Substitution of this formula into the conservation laws (2.1), (2.2) yields the equations Finally, we show that the flow rule (2.18) may be rewritten To see this, we combine (2.24b) with (2.25) to conclude and then substitute the relation α(I) = τ/p + p/C(φ) derived by manipulating ( . (2.37) Substituting (2.37) into the continuity equation (2.35a), we find If φ = 0, then this equation reduces to div u = 0, although this is of course a highly singular limit. Thus, if φ = 0, the left-hand side of (2.36) vanishes, so this equation simplifies to the yield condition τ = μ(I)p, and substitution into (2.35b) yields (2.15). This proves the lemma.

Proofs, part I: linearization (a) An alternative formulation of the alignment condition
It is convenient to study the linearized equations with a reformulated alignment condition that describes stress in terms of eigenvectors of, rather than entries of, the stress tensor. Since τ defined by (2.4) has trace zero, it has eigenvalues 6 ± τ . Taking ψ as the angle that the eigenvector with which may be verified by checking that (cos ψ, sinψ) is an eigenvector of this matrix with eigenvalue − τ . Thus, the stress tensor σ ij is completely specified by the three scalars p, τ and ψ.
Focusing on the first rows of the strain-rate tensor (2.9) and of (3.1), we extract from the matrix equation (2.8) the vector equation where k = −2 D < 0. Since D and τ lie in the two-dimensional space of trace-free, symmetric matrices, (3.2) is equivalent to (2.8). It follows from (3.2) that Alt. alignment: In point of fact, this equation is slightly weaker than the alignment condition since (3.3) is consistent with the possibility that k > 0 in (3.2); to rule out the latter possibility we impose the supplemental inequality 7 that Substitution of the stress tensor (3.1) into the momentum balance equations (2.2) allows for the full set of equations to be written as and This system has five scalar unknowns, U = (u 1 , u 2 , φ, p, ψ). In (3.5a), (3.5b), τ is a mnemonically suggestive abbreviation for the yield function Y(p, φ, I) in (2.17), and in (3.5d), a repetition of (2.18), the function f depends on arguments (p, φ, I) that are not written explicitly. As in appendix B, to linearize the equations we substitute a perturbation of a base solution into the equations, retain only terms that are linear in the perturbationÛ and freeze the coefficients at an arbitrary point (x * , t * ). It is convenient to temporarily drop most terms not of maximal order and estimate their effect in a calculation at the end of the argument. For example, this construction applied to (3.5c) yields the constant-coefficient, linear equation . Lower-order terms ∂ j φ * û j and ∂ j u * jφ in the full linearization of (3.5c) have been dropped in (3.7).
In expanding the fully nonlinear factor D in (3.5d), we may take advantage of the rotational invariance of the equations to arrange that ψ * = 0; i.e. we may calculate in a rotated coordinate , the x 1 -axis is the maximal stress axis. Then by the alignment condition (3.3), the base-state deviatoric strain-rate tensor is diagonal at (x * , t * ) and by (3.4), in the 1,1-position of this matrix, ∂ 1 u * 1 − ∂ 2 u * 2 < 0. This corresponds to non-zero compression along the major stress axis, as illustrated in figure 3. Now where the approximation follows from the expansion if A > 0 and |X|, |Y| A. Thus, as given in table 2, the (local) linearization of D equals −(∂ 1û1 − ∂ 2û2 )/2. In (3.5d), the function f contains p, φ and I as implicit arguments. As reflected in the table, the dependence on p and φ contributes zeroth-order terms in these variables to the linearization.
In (3.5a), (3.5b), τ also depends on p, φ and I, and the terms involving τ are differentiated; hence new issues arise in linearizing them. For example, by the chain rule,    Since ψ * = 0, the full linearization of, say, the first term here equals (∂ p τ ) * ∂ jp , a term given in the table, plus lower-order terms All of these terms, as well as numerous other analogous terms in the full linearization of (3.5a) that are not of maximal order, have been dropped in (3.10a)-(3.10e). Putting all the pieces together, we obtain the linearization 8 of the system (3.5a)-(3.5e) and and Observe that, by hypothesis (2.19), B = C, a fact that we use in (4.3) and below. 8 These equations are maximal order except that in (3.10a) and (3.10b) the term d * tûj retains first-order spatial derivatives even though these equations also contain second-order derivatives ofû j . whereŨ = (ũ 1 ,ũ 2 ,φ,p,ψ) is a 5-vector of scalars, ξ = (ξ 1 , ξ 2 ) is a vector wavenumber, , indicates the inner product and λ is the growth rate. The function (4.1) is a solution of (3.10a)-(3.10e) iff λ,Ũ satisfies the generalized eigenvalue problem where u * = (u * 1 , u * 2 ), On the right side of (4.2), the modified eigenvalue parameter is λ + i u * , ξ because where and S 12 , S 21 and S 22 fill out the rest of the matrix. DefiningŨ 1 = (ũ 1 ,ũ 2 ,φ) andŨ 2 = (p,ψ), we rewrite (4.2) as The zero entries in the last two rows of E mean that S 21Ũ1 + S 22Ũ2 = 0 so we can solve for Substitution ofŨ 2 into (4.7) then reduces this problem 9 to the ordinary 3 × 3 eigenvalue problem, where E 11 is the 3 × 3 block in the upper left of E. 9 In other words, we are performing on the symbol level the reduction that we performed on the operator level in §2d. We decompose the 3 × 3 matrix in (4.9) into smaller blocks, ⎡ where we calculate as the contribution of S 11 , as the contribution of −S 12 S −1 22 S 21 , which is symmetric, and (4.13)

(b) Estimation of the eigenvalues
We claim that the growth-rate eigenvalues (4.10) satisfy Since only the real parts of eigenvalues matter, we may drop the term i u * , ξ in (4.10) and verify (4.14) for the eigenvalue problem 10 where we write for the matrix in (4.10) and we shorten the notation by dropping the subscript 1 onŨ. For large ξ , it is instructive to use perturbation theory to compare the eigenvalues (4.15) with the eigenvalues P 0Ũ = −ΛŨ, where accounts for the cross terms. For the specific matrices (4.11) and (4.12), det M = 0, and This proves the lemma.

Remark 4.2.
It is noteworthy that det N > 0 except for the two directions Effectively, this calculation rederives the result of Pitman & Schaeffer [12] that the equations of CSSM, even without I-dependence, are well-posed for all directions except possibly those defined by (4.22).
It follows from lemma 4.1 that P 0Ũ = −ΛŨ has two eigenvalues, say Λ 1 , Λ 2 , where Λ 1 , Λ 2 < 0 and is homogeneous of degree 2 in ξ . Since P is an O(|ξ |)-perturbation of P 0 , two of the growthrate eigenvalues of (4.15) satisfy both of which are negative in the limit |ξ | → ∞; i.e. they are bounded above by zero in this limit. The third growth rate is given by The first term on the extreme right is the ratio of two quartics, the denominator being non-zero, so it is bounded, and the perturbation decays at infinity. This verifies (4.14) for all three eigenvalues derived from (3.10a) to (3.10e). It remains to consider the effect of the lower-order terms that were neglected in (3.10a)-(3.10e). Inclusion of these terms would lead, after a calculation as above, to an eigenvalue problem (4.15) for a perturbed matrix ⎡ As above, two of the eigenvalues of this matrix are negative and O(|ξ | 2 ), and invoking the determinant shows that the third is bounded. This verifies (4.14) for eigenvalues of the full linearization of (3.5a)-(3.5e) and hence shows that the system is linearly well-posed. It is, therefore, expected that numerical solutions of the full two-dimensional nonlinear transient equations will not exhibit the exponential blow-up of perturbations or the dependence on grid resolution seen in the incompressible μ(I) equations [9].

Chute flow in compressible I-dependent rheology
Let us recall the steady-uniform 'Bagnold' solution [4,31] to the (incompressible) μ(I)-rheology [4,5] for chute flow. Assuming that μ(I) is given by (2.13), these solutions exist for inclination angles ζ between the maximum and minimum angles ζ 2 = tan −1 (μ 2 ) and ζ 1 = tan −1 (μ 1 ), respectively. Letting Oxz be a Cartesian coordinate system, with the x-axis pointing downslope and the z-axis being the upward pointing normal, the Bagnold solution is where h is the flow depth, φ is a constant solids volume fraction and u is the downslope velocity. At a fixed inclination, the inertial number is equal to the positive constant There is strong experimental and discrete element method (DEM) evidence [4,5,31,32] for both the lithostatic pressure distribution and the three-halves power in the dependence of the velocity on the thickness in (5.1). Fortunately, the CIDR solution is close to this, as we now show. Assume that all variables depend only on z and only the x-component of the velocity is nonzero, say u = u(z). Motivated by Bagnold flow, the scalings are used to non-dimensionalize the variables and the following ordinary differential equations (ODEs) are derived forp andǔ below.

Lemma 5.1. Under the above assumptions, the non-dimensional pressurep(z) and the non-dimensional downslope velocityǔ(z) predicted by the CIDR rheology satisfy
and dǔ dž = p, (5.5) where the pressure is zero at the free surfacep(1) = 0, the velocity is zero at the baseǔ(0) = 0 and for the function C(φ) proposed in (2.30) and (2.31) the solids fraction φ = φ(p) is Steady uniform solutions to the CIDR model are shown in figure 4 for χ = 0.5, 2, 5, 20, 100 and with φ max = 0.6 and φ min = 0.5. The nature of the solution is controlled by the non-dimensional parameter χ , which is larger for thicker flows. This parameter is also inversely proportional to β(I ζ ), which by (5.2) and (2.27) is a weakly increasing function of the inclination angle as shown in figure 5. For large χ , the concentration graphed in figure 4a is close to φ max over a significant proportion of the flow depth, with a narrow boundary layer near the free surface, where the concentration decreases to φ min . In the limit as χ → ∞, the solution tends to the red line in figure 4a, which corresponds to Bagnold flow with φ = φ max . As χ is decreased, the surface boundary layer becomes thicker and the flow becomes progressively more dilute, tending to the blue line in figure 4a as χ → 0, which corresponds to Bagnold flow with φ = φ min . In figure 4b,c, the non-dimensional pressure and the velocity are shown, which are also bounded between the red and blue lines corresponding to Bagnold flow with φ equal to φ max and φ min , respectively. As well as the profiles being bounded by two non-dimensional Bagnold solutions, the steadyuniform CIDR solutions also have exactly the same scaling properties (5.3), on the flow density ρ , gravity g, the flow depth h, the chute inclination ζ , the particle diameter d and the inertial   and dp dz + ρ * φg cos ζ = 0. (5.8) We may apply the chain rule to (5.7) and (5.8) to derive dτ/dp = tan ζ , and, since tractions at the surface vanish, τ = p tan ζ .

Conclusion and discussion
In this paper, we have analysed a generalization of the μ(I)-rheology that allows for changes in the granular solids volume fraction. The equations of motion, the CIDR model, are found to be linearly well-posed when the constitutive laws satisfy certain criteria. We indicate (in §2e) how the specific I-dependence of the μ(I)-rheology can be fitted into the theory, and we have shown that the inclusion of compressibility removes the ill-posedness of the incompressible rheology [9]. CIDR equations for steady, uniform flow down an inclined chute are solved in §5, where we observe that the solution is comparable to the classical Bagnold solution under conditions similar to the many experimental and DEM results [4,5,31,32] for the velocity profile in steady chute flows. Crucially, the equations are well-posed at all inertial numbers, which corresponds to the full range of inclination angles for steady chute flow, in contrast with the limited range of angles in which the incompressible theory is well-posed [9]. At the same time, the new theory captures the dilatant behaviour of over-consolidated granular material as the structure of the equations is motivated by CSSM.
Incidentally, the chute flow calculation relates to Forterre & Pouliquen [33], who modify incompressible μ(I)-rheology by postulating a constitutive law in which φ is explicitly specified as a function of the inertial number, say Φ(I). In the CIDR model, there is no constitutive law relating these two variables; nevertheless, dependence of φ on I is implicit in the solution. It is also quite possible that a more subtle dependence of volume fraction on inertial number will be needed to reproduce some phenomena reported recently in experiments and numerical simulations, while retaining the property of well-posedness.
Our primary viewpoint in this paper has been to regard the CIDR model as modifying μ(I)rheology with compressibility. However, it is equally valid to regard CIDR as modifying CSSM with rate dependence. To recapitulate, our result shows well-posed equations result from such modification provided that the yield locus and flow rule satisfy (2.19) and (2.20). Unlike in the incompressible μ(I)-rheology and the rate-independent CSSM equations, the CIDR equations are linearly well-posed for all deformations and for perturbations in all directions in Fourier space.
An important next step will be to specify constitutive laws, satisfying the general conditions for linear well-posedness of this paper, and formulated to accurately match the available experimental results and discrete numerical simulations for granular flows such as twodimensional steady flow in a Couette geometry and time-dependent chute flow. Then the CIDR model and well-posedness result can be tested with fully two-dimensional nonlinear transient numerical computations of these flows. We refer to Jackson [1] for a derivation of (A 2) from the normality condition of plasticity. By way of example, a simple, physically acceptable, yield locus is given by where μ is a coefficient of friction, as in (2.5), and C(φ) is an increasing function of the solids fraction of the form (2.30). For such a yield condition, it follows from the proof of lemma 2.5, restricted to the case where μ(I) is independent of I, that equations (2.1), (2.2), (A 1), (A 2), (2.8) reduce to the Coulomb model in the limit φ → 0.

(b) Consequences of the flow rule
The behaviour discussed in this subsection occurs under fairly general hypotheses-see Jackson [1]. However, to explain the theory with a minimum of technicalities, we confine the discussion to the specific yield condition (A 3). The phrase critical state, from which CSSM derives its name, refers to a state p, φ such that According to the flow rule (A 2), at a critical state, deformation is not accompanied by any change in φ. Let us examine behaviour away from the critical state line. Suppose that, for example, initially the (uniform) state of material is at yield at point A in figure 2b. At this point, ∂Y/∂p < 0, so according to the flow rule div u < 0; i.e. material compactifies and becomes stronger, so τ must increase for deformation to continue. Indeed, the stress will continue to increase until a critical state on a larger yield surface is reached, as suggested in figure 2b by the φ 3 -yield surface. Moreover, if φ in (2.30) is small, a very slight increase in φ is sufficient to accommodate this evolution. That is, we expect stress to be quickly driven from point A to a critical state on a larger yield surface where the Coulomb yield condition (A 6) is satisfied.
Conversely, at point B in figure 2b, ∂Y/∂p > 0, so under deformation div u > 0; i.e. material expands and becomes weaker. It is natural to imagine that the stress is driven to a critical state on a smaller yield surface, as suggested by the arrow in the figure. This would indeed be the case if material deformed uniformly, but this assumption is unrealistic for stresses above the critical state line, τ > μp. For such stresses, because material expands under deformation and therefore weakens, instability often causes localized deformation-if deformation near one point happens to be slightly larger than elsewhere, the associated expansion lowers the yield condition more near this point, and subsequent deformation tends to concentrate near this point. We illustrate the above test on the following made-up nonlinear system that has some similarity to the PDEs analysed in this paper: The linearized equations with frozen coefficients are where u * , v * is the base-state solution evaluated at the point (x * , t * ) andû,v are the perturbations solutionsÛ whereŨ ∈ R 2 satisfies the eigenvalue problem The eigenvalues of (B 7) could be easily calculated exactly but, provided that u * = 0, they can be estimated more easily from their asymptotic behaviour as |ξ | → ∞, Note that, in analysing linear ill-posedness of (B 4), we consider the full linearization of the equations, i.e. (B 5). One might be tempted to discard terms with lower-order derivatives in the expectation that the growth of exponential solutions as |ξ | → ∞ ought to be dominated by the highest-order derivatives in the equation. However, the counter-example shows that this expectation is not valid, in general.
Nevertheless, for this example we may in fact analyse exponential solutions of (B 5) by first considering the maximal-order linearized equations In each of the above equations, only terms of maximal order are retained, i.e. the terms (∂ xx v) * û andv have been dropped from the second equation because it contains the higher-order terms ∂ xû and u * ∂ xxv , respectively. The growth rate of exponential solutions of (B 10) satisfy the same estimates (B 8), and the neglected lower-order terms do not change the leading-order behaviour. Often calculations may be simplified by studying the maximal-order linearization as an intermediate step.

(b) Consequences of ill-posedness (i) Restrictions on the existence of solutions
In order for an ill-posed initial value problem to have a solution, usually the initial conditions must satisfy an extreme smoothness requirement, stronger than is physically acceptable in most applications. It may be difficult to demonstrate this behaviour, in general, but, for the backwards heat equation, we illustrate the behaviour with the following.
Proposition B.1. If ε = 0, the initial value problem for (B1) with initial data u(x, 0) = ε| sin x| p has no solution for any positive time interval unless the power p is an even non-negative integer. 12 Proof. Suppose the 2π -periodic function f (x) has Fourier series f (x) ∼ c n e inx . If equation (B 1) with initial condition u(x, 0) = f (x) has a continuous solution for 0 ≤ t < η, it has the Fourier series representation Now if the Fourier coefficients of a function f (x) satisfy n 2k |c n | 2 < ∞, then the derivatives (d/dx) j f (x) are square integrable for j = 0, 1, . . . , k. But, provided p is not an even non-negative integer, the proposed initial data | sin x| p has singular behaviour near x = 0. Specifically, (d/dx) k | sin x| p is square integrable only if k < p + 1/2. It follows for the Fourier coefficients of | sin x| p that, if k > p + 1/2, This inequality is incompatible with (B 11), so the initial value problem cannot be solved on any positive time interval.

(ii) Grid-dependent computations
The attempt to solve an ill-posed PDE numerically produces unreliable, grid-dependent, results. Such behaviour has been observed in various physical problems [7][8][9][10] where the formulation was based on an ill-posed system of equations. However, in complicated problems like these, usually computational resources are stretched to the limit, meaning behaviour under grid refinement cannot be readily probed. Let us illustrate grid dependence on a much less demanding problem, the toy problem (B 4) above.
If ε = 0 and with initial conditions  is also a steady-state solution of the nonlinear system (with ε > 0). If a > 1, then u(x, ∞) > 0, and the calculations above suggest that the equations will be linearly well-posed. However, if a < 1, then u(x, ∞) dips below zero over an interval, which suggests that the equations will be ill-posed. Figure 6, where a = 1.5, and figure 7, where a = 0.5, confirm these expectations. They show numerical solutions using a central-space forward-time explicit scheme on the periodic domain x ∈ [−π , π ] with a spatial resolution of x = 2π/100. Each figure has two panels, one showing the temporal evolution of the distance from the asymptotic solution, ) and the other plotting the two variables u, v at a specific (late) time during the computation. In figure 6, the well-posed case, the numerical solution converges to the predicted steadystate solution (B 14) within numerical accuracy. By contrast, in figure 7, after an initial decay, ill-posedness asserts itself and causes the solution to blow up. Regarding grid dependence, figure 8 shows another computation in the ill-posed case with a coarser grid, x = 2π/30. The solution appears to converge to the steady-state solution, just like in the well-posed case. In other words, the computations on the coarse grid hide the ill-posed character of the underlying PDEs. This highlights that, in order to extract meaningful information from numerical computations, a proper study of grid convergence must be first carried out. Incidentally, if the grid is made finer than in figures 6 and 7, in the well-posed case a > 1 the numerical solution converges to the steady-state solution with a smaller numerical error, while in the ill-posed case it blows up sooner, as expected.