3.1. Uniform geometry

Figure 2 shows the dimensionless flux per unit length, $Q$ , as a function of the dimensionless pressure drop, $\Pi$ , for planar (i.e. independent of $x$ ) steady-state solutions in the uniform geometry, $h=1$ , calculated using pseudo-arclength continuation. Here, we vary the value of the parameter governing the viscosity contrast, $\gamma =\log (\mu _c/\mu _h)$ . In the following discussion, we report the quantitative results of the current model system, while a direct comparison to the corresponding values obtained from the models of Wylie & Lister (Reference Wylie and Lister1995), Helfrich (Reference Helfrich1995) and Morris (Reference Morris1996) is given in Appendix B. The qualitative behaviour of the system is identical to that of the models in these previous studies. Namely, at low values of the viscosity contrast, the solution curve is single-valued and the planar steady states are stable. Above a first critical value, $\gamma _c=4.8$ , the curve becomes multivalued, resulting in three branches of the solution curve: a ‘hot’, fast branch at high pressure drops; a ‘cold’, slow branch at low pressure drops; and an intermediate branch on which the flux increases with decreasing pressure drop. As in previous studies, the intermediate branch is found to be unstable to a planar instability (i.e. to perturbations with vanishing wavenumber in the $x$ -direction), which acts to drive solutions off the intermediate branch and towards the fast or slow branches of the solution curve. Above a second, higher, viscosity ratio, $\gamma _{3d}\approx 5.2$ , there exist solutions which are most unstable to perturbations with non-vanishing wavenumber in the $x$ -direction (see the dotted region of the solution curves for $\gamma =5.37$ and $\gamma =5.5$ in figure 2). This thus marks the onset of a fingering instability, analogous to classical Saffmann–Taylor viscous fingering, by which the uniform steady states go unstable to the formation of alternating regions of hot, low-viscosity fluid, and cold, high-viscosity fluid. At $\gamma$ just above $\gamma _{3d}$ , the region of the solution curve which is most unstable to the fingering instability is limited to a small region near the fold bifurcation that joins the slow and intermediate solution branches. As $\gamma$ is increased further, this region extends down the slow solution branch and the entire slow branch becomes unstable to the fingering instability at a third critical value of the viscosity ratio, $\gamma _\infty \approx 5.4$ . Thus, by $\gamma =5.5$ , the entire slow branch is unstable to the fingering instability and, as we will discuss later, we are able to find stable non-planar steady states reached after the onset of fingering. We report approximate values of $\gamma _{3d}$ and $\gamma _\infty$ here, because these results were calculated using our numerical method, with a finite domain ( $L_x=1$ in the cases shown in figure 2) and symmetry boundary conditions. This does not affect the calculation of $Q$ for the planar steady states, nor the stability of the solutions to the planar instability; however, it does impact the stability to non-planar perturbations, since the finite domain restricts perturbations to those fitting inside the domain (with a half-integer number of wavelengths). Thus, for an infinite domain, $\gamma _{3d}$ is likely slightly smaller than reported here. Similarly, the finite resolution of the grid sets a lower bound on the wavelengths that can be accurately captured. Since the wavelength of the most unstable mode decreases with decreasing $\Pi$ (see later and Wylie & Lister Reference Wylie and Lister1995; Morris Reference Morris1996), determining $\gamma _\infty$ (where the fingering instability first extends to arbitrarily small values of $\Pi$ ) is difficult with our numerical method. Ultimately, our aim is to study the impact of a spatially varying channel width at a fixed value of $\gamma$ and moderate values of $\Pi$ (in the neighbourhood of the fold bifurcation), and not to accurately capture these critical values of $\gamma$ for an infinite uniform geometry.

Figure 2.

Dimensionless flux per unit length, $Q$ , against dimensionless pressure drop, $\Pi \equiv \Delta p h_0^4/\kappa \mu _h L^2$ , for planar steady states in the uniform geometry, $h=1$ . Results are shown for several values of $\gamma =\log (\mu _c/\mu _h)$ . The dimensional flux per metre is given by $\kappa L Q/h_0$ . Solid lines show stable regions of the solution curve, and dashed and dotted lines show regions that are most unstable to planar and non-planar perturbations, respectively. See table 2 for other parameters used in these solutions.

We now provide some more detail of the steady-state solutions for the uniform geometry, $h=1$ , and at $\gamma =5.5$ , fixed (now on a larger domain, $L_x=10$ ). Figure 3(a) again shows the total flux per unit length of the fissure, $Q$ , as a function of $\Pi$ at steady state. The black (solid, dashed and dotted) line shows the solution branch for the planar solutions (independent of $x$ ), as in figure 2. Point A indicates the fold bifurcation, below which the planar steady states are unstable. The wavenumber, $k$ , of the most unstable mode for the unique planar steady state with flux, $Q$ , is shown in figure 3(b). In panels (a) and (b), point B indicates the location at which the most unstable mode becomes non-planar ( $k\neq 0$ ). Assuming the eruption starts at a high pressure drop and wanes with time, the structure of the bifurcation diagram gives the potential for the system to pass through the fold bifurcation and exhibit a sudden reduction in flux before localising due to the fingering instability. As reported by Wylie & Lister (Reference Wylie and Lister1995) and Morris (Reference Morris1996), the wavenumber of the most unstable mode increases as $\Pi$ decreases, and scales with $1/Q$ (see the red dashed curve in figure 3 b) since the length scale is set by the thermal entry length. For our results, the wavenumber curve in figure 3(b) is stepped, which is simply the signature of the finite domain (in this case, $L_x=10$ ) and symmetry boundary conditions, meaning that the wavenumers can only take half-integer multiples of $1/L_x$ – however, by taking the domain size to infinity, we could retrieve a perfectly continuous spectrum of wavenumbers.

Figure 3.

(a) Dimensionless flux per unit length of the fissure, $Q$ , as a function of dimensionless pressure drop, $\Pi$ , for stable steady states with $h=1$ and $\gamma =5.5$ (as in figure 2). The colour indicates the surface flux ratio, $Q_r$ defined by (2.29), and the planar steady-state solution curve is shown by the black dashed and dotted line. Point A shows the fold bifurcation at which the solutions go unstable to planar perturbations, and point B indicates the location at which the solutions become most unstable to non-planar perturbations. (b) Wavenumber of the most unstable mode, $k$ , as a function of $Q$ (black). The red dashed line shows a curve proportional to $1/Q$ . Points A and B as in panel (a). (c) Dimensionless cross-fissure average temperature, ${\mathcal{T}}$ defined in (2.13), for the last stable steady state near the slow branch at $\Pi =264$ . The domain has been reduced to show a single wavelength of the fingering pattern (of wavelength $0.909$ ). The contours are plotted from ${\mathcal{T}}=0.1$ near the top and increase in increments of 0.1. (d) Contour plot of vertical flux, $q_z=2\Pi {\mathcal{L}}\partial p/\partial z$ , for the solution shown in panel (c). The central (closed) contour is for $q_z=2.4$ , and values decrease in increments of 0.2 moving out from here. Note that the results in $x\lt 0$ have been produced by reflecting in $x=0$ and have not been explicitly calculated. See table 2 for other parameters used in these solutions.

Figure 3(a) also shows coloured points representing stable steady states. As noted previously, the fast branch of the planar solutions is stable. The stable solutions near the slow branch were obtained by natural continuation, starting at $\Pi =200$ and increasing $\Pi$ in steps of 0.5. At each step, the new steady state is either stable, in which case, it is kept; or else, it is unstable, in which case, it is perturbed by its most unstable perturbation and evolved forwards in time until a new steady state is reached. The colour is set by the surface flux ratio, $Q_r$ , defined by (2.29). Thus, $Q_r=1$ on the fast, planar branch, and $Q_r\gt 1$ near the slow branch where there are flow-focussed regions. We found that at values of $\Pi$ close to the onset of the non-planar instability (point B in figure 3), the perturbed solutions evolve onto the fast branch of the solution curve, resulting in no non-planar steady states, as was observed by Wylie & Lister (Reference Wylie and Lister1995) when trying to find non-planar steady states. However, at lower $\Pi$ , we find that the non-planar solutions evolve towards stable steady states with localised finger regions, as suggested by the solutions of Helfrich (Reference Helfrich1995). This result was also reproduced for the unaveraged system of Wylie & Lister (Reference Wylie and Lister1995), but at a higher viscosity ratio, $\gamma =6\gt \gamma _\infty$ , such that the entire slow branch is unstable to the fingering instability. The total flux of these flow-focussed states is slightly increased compared with the unstable steady state, and is distributed into higher and lower flux regions along the fissure. Among these solutions, the surface flux ratio and the wavelength of the finite amplitude fingering pattern both increase with $\Pi$ , reaching a maximum surface flux ratio of $Q_r\approx 1.8$ at the last pressure drop before the non-planar solutions become unstable and instead evolve onto the fast branch (in this case, $\Pi =264$ ). We again note that on this steady flow-focussed branch, the finite domain size and symmetry boundary conditions are relevant, since only a half-integer number of wavelengths can fit in the domain. Thus, the increase of the wavelength with increasing $\Pi$ does not occur continuously, but in discrete jumps marked by a sequence of bifurcations. Indeed, since the solutions are calculated by natural continuation, figure 3(a) does not capture all stable steady-state branches; rather, there is a separate solution branch corresponding to each half-integer number of fingers within the domain, and the regions over which each branch is stable are overlapping. Thus, in figure 3(a), we follow a branch up to the $\Pi$ at which it goes unstable, and then land on the next stable branch which is again captured up until the $\Pi$ at which it goes unstable, and so on. The bifurcation which results in the onset of instability on one branch, and the transition to the next, is a subcritical pitchfork bifurcation (e.g. see Strogatz Reference Strogatz2015a ), involving the symmetry-breaking selection of which end of the domain switches from a hot to a cold finger (or vice versa). In contrast, for an infinite domain, $L_x=\infty$ , we would anticipate that the flow-focussed steady-state branch would be a single continuous curve with a continually varying wavelength.

Figure 3(c,d) shows contour plots of the cross-channel average temperature, ${\mathcal{T}}$ , and vertical flux, $q_z$ , for the final stable steady solution on the non-planar steady-state curve at $\Pi =264$ , evidencing a finger structure with a wavenumber of $k=1.1$ (i.e. 11 wavelengths in the domain of width $L_x=10$ ), which is close to the most unstable wavenumber, $k=1.2$ , predicted for the uniform steady state at the same pressure drop. This figure also demonstrates why the ratio of fluxes at the surfaces, $Q_r=1.8$ , is not particularly large (compared with values obtained for the non-uniform geometry discussed later, see § 3.2), since the localisation occurs most significantly at depth before becoming more uniform nearer the surface where the magma is close to the wall temperature, $T=0$ . Under the model of a waning eruptive pressure drop (and assuming that the pressure drop varies sufficiently gradually that the system progresses through steady states), the system would drop off the fast branch at the fold bifurcaton around $\Pi =240$ and first attain a non-planar solution at this value. In this case, the surface flux ratio would be smaller, resulting in a $14\,\%$ difference in flux between the hot and cold regions of the outflow. It should be noted that at $\gamma =5.5$ , we are not particularly far above the onset of the fingering instability (at $\gamma _{3d}\approx 5.2$ ), and one would anticipate that at larger viscosity ratios, the steady-state fingered solutions would exhibit a greater degree of flow-focussing. This will also be true of the geometrically flow-focussed states, discussed in the following sections, and which we will show are associated with the fold instability arising first at $\gamma _c=4.8$ .

In this section, we have revisited the uniform geometry case, and shown that our results provide close qualitative agreement with the previous findings of Helfrich (Reference Helfrich1995), Wylie & Lister (Reference Wylie and Lister1995) and Morris (Reference Morris1996), evidencing multiple steady states and the onset of the thermoviscous fingering instability at low pressure drops. Crucially, these results provide the baseline case against which to interpret the results in the non-uniform geometry, which we consider next. A quantitative comparison to the alternative averaging approaches of Helfrich (Reference Helfrich1995), Wylie & Lister (Reference Wylie and Lister1995) and Morris (Reference Morris1996) is made in Appendix B, showing that the heat balance averaging of § 2.1 improves on the averaging of Helfrich (Reference Helfrich1995), but that there are nonetheless quantitative discrepancies with the unaveraged system. We note that there are multiple other approximations that are likely to have equally significant effects on the quantitative predictions of all models, in particular, the simple treatment of the thermodynamics at the fissure walls and the particular choice of the viscosity function, (2.5). For example, results of Wylie & Lister (Reference Wylie and Lister1995, p. 248) employing a step-function model for the viscosity find a critical viscosity ratio, $\mu _c/\mu _h$ , an order of magnitude smaller than when employing the exponential viscosity model (12 compared with 190), and Morris (Reference Morris1996, p. 126) finds a critical viscosity ratio of just 7.3 when the temperature of the fissure walls vary linearly with flow direction. In using the heat balance averaging detailed in § 2.1, we thus accept an additional source of approximation and associated quantitative discrepancy with the unaveraged system, in exchange for a number of advantages. As noted before, in employing the heat balance averaging, we consistently average the time derivative in the heat equation, meaning that time-dependent states (for example, those discussed in §§ 3.4 and 4) can be meaningfully considered alongside the steady states of the system – with both arising from the same assumed structure of the temperature profile and resulting average of the heat equation. Second, from a purely practical standpoint, in the unaveraged approach, evaluation of the functions ${\mathcal{L}}$ and $\mathcal{S}$ , and their derivatives, relies on interpolation and numerical differentiation of a precomputed table of values, $\bar {\mu }(\tau )$ . Furthermore, the derivative $\textrm {d}{\mathcal{L}}/\textrm {d}\tau$ diverges at $\tau =0$ . These features make the numerical solution of the resulting equations harder than when employing the analytical expressions for ${\mathcal{L}}$ and $\mathcal{S}$ arising from the heat balance averaging.

Figure 4.

Steady-state solution curves for different amplitudes and wavelengths of the dimensionless fissure half-width, $h=1+A\cos (2\pi x/\lambda )$ . Each panel shows average flux per unit length, $Q$ , as a function of the dimensionless pressure drop, $\Pi$ , as in figure 2, with points coloured by the scaled surface flux ratio, $\tilde {Q}_r$ , defined by (2.30) (note the logarithmic colour scale). The wavelength, $\lambda$ , is shown in the title and the three coloured curves correspond to amplitudes of $A=0.02$ , $0.05$ and $0.1$ (increasing in the direction shown). The solution curve for the planar steady states in the uniform geometry ( $A=0$ ) is shown as the black dotted curve in all panels. See table 2 for other parameters used in these solutions.

Figure 5.

Colour plot of scaled surface flux ratio, $\tilde {Q}_r$  (2.30), as a function of dimensionless pressure drop, $\Pi$ , and the wavelength of the variations in width, $\lambda$ , at different values of the amplitude, $A$ . Contours (black) are plotted at flux ratios of $\tilde {Q}_r=4$ and 8 in all panels, and at the additional level $\tilde {Q}_r=16$ , in panel (d). The axis ticks indicate the location of solutions used to construct the contour plot, and red dots indicate that an average over the period of a periodic solution was used as opposed to a steady-state solution. Note the extended $\Pi$ -axis in the final panel. See table 2 for other parameters used in these solutions.

3.2. Non-uniform geometry

We now introduce a sinusoidal perturbation to the fissure width, at different amplitudes and wavelengths. Figure 4 shows the corresponding $Q\! (\Pi )$ steady-state solution curves for $\lambda =0.4, 1$ and $2$ , and $A=0.02, 0.05$ and $0.1$ . The solution curve corresponding to the planar solutions of the uniform geometry ( $A=0$ ) is shown as a dashed line and the colours now indicate the scaled surface flux ratio, $\tilde {Q}_r$ defined by (2.30). We also now use a logarithmic scale for the colour bar. We see that the effect of the non-uniform geometry is generally to replace the unstable intermediate branch of the planar solutions with a stable solution branch, exhibiting significant flux focussing. This stable, ‘focussed’, branch is either part of a continuous single-valued solution curve, or is connected to the original fast and/or slow branches by one or two pairs of fold bifurcations (for example, for $A=0.05$ and $\lambda =2$ , there are two pairs of fold bifurcations, one at either end of the stable central branch). The deviation from the uniform geometry solution increases with both amplitude and wavelength. Further results for the scaled surface flux ratio, $\tilde {Q}_r$ , are shown for a wider range of wavelengths and amplitudes in figure 5. For this plot, the solutions were found by natural continuation in the pressure drop, $\Pi$ , starting on the fast branch with $\Pi =400$ and integrating forwards in time to obtain stable steady states as $\Pi$ is decreased in discrete steps. This clearly evidences the existence of the highly localised solutions at intermediate pressure drops, and the dependence on the wavelength and amplitude of the fissure geometry. There are a small number of parameter values, primarily lying on the slow solution branch, at which no steady state was found (indicated by red dots in figure 5), and the solutions were instead found to be periodic. Oscillatory dynamics have previously been observed by Pailha et al. (Reference Pailha, Hazel, Glendinning and Juel2012) in experiments on isothermal fingers of air displacing viscous fluids, where geometrical constraints were found to play a key role (see also Thompson, Juel & Hazel Reference Thompson, Juel and Hazel2014; Franco-Gómez et al. Reference Franco-Gómez, Thompson, Hazel and Juel2016; Lawless, Hazel & Juel Reference Lawless, Hazel and Juel2024). The structure of the periodic solutions observed in these cases is rather different from those observed here. We explore the periodic solutions in more detail in § 3.4, but note that since these have a period-averaged surface flux ratio in line with neighbouring steady solutions (see figure 5), the implication for the degree of fissure localisation is minimal. The temperature and flux distributions are shown in figure 6 for steady-state solutions at $\Pi =240$ on the focussed branch for $A=0.05$ and $\lambda =2$ , 1 and 0.5. These can be compared with figure 3(d,e) for the focussed solutions arising from thermoviscous fingering in a uniform geometry at $\Pi =264$ , showing that the temperature variations penetrate further up the fissure and that the vertical flux variations are several times larger.

Figure 6.

(a) Cross-channel average temperature, ${\mathcal{T}}$ , defined in (2.13), for steady-state solutions at $\Pi =240$ with $A=0.05$ and $\lambda =2,\,1,\,0.4$ (from left to right). The colour bar is shared between panels, and the contours are plotted from ${\mathcal{T}}=0.1$ near the top left/right and increase in increments of 0.1. (b) Vertical flux, $q_z\equiv 2{\mathcal{L}}\partial p/\partial z$ , for the same solutions. Contours are plotted from $q_z=2$ at the left and right, and increase inwards in increments of 2. Note that the results in $x\lt 0$ have been produced by reflecting in $x=0$ , and have not been explicitly calculated. See table 2 for other parameters used in these solutions.

The physical origin of the stable, highly localised branch can be readily understood as the wider region of the fissure remaining in the hot, fast flowing state, while the thinner region transitions to the cold, slow state as the pressure drop decreases. The flow in the wider region is cognate with the fast branch of the planar solution, and the flow in the narrower region is cognate with the slow branch of the planar solution. Thus, the occurrence of geometrical focussing is more closely associated with the fold bifurcation and planar instability exhibited by the uniform-width system than it is with the thermoviscous fingering instability. The dependence on amplitude and wavelength can also be easily explained. Naturally, a larger amplitude results in a greater difference in flow rate and hence heat advection through the different regions of the fissure. However, if the wavelength is small, then the hot and cold regions can exchange heat more easily, resulting in a more uniform temperature and a reduction in the geometric effect. Note that this exchange of heat does not occur diffusively, since the reduced heat equation includes no diffusion in the plane of the fissure (see (2.9)); rather, it occurs through advection due to along-fissure pressure gradients. These pressure gradients are stronger when the length scale between the regions of the fissure are smaller, and thus the solutions become more uniform for smaller $\lambda$ (e.g. compare solutions for $\lambda = 2$ and $\lambda = 0.4$ in figure 6). Another way to understand the suppression of geometrical localisation at short geometrical wavelength is by reference to the fingering instability in the uniform geometry. Here, Wylie & Lister (Reference Wylie and Lister1995) and Morris (Reference Morris1996) showed that the preferred wavelength of the fingering instability is comparable to, or larger than, the thermal entry length, which scales with $Q$ . The geometrical localisation occurs via the above mentioned mechanism at values of $\Pi$ near the location of the fold bifurcations associated with multiplicity of steady states in the planar problem (see figure 4), and hence at solutions which have fluxes of the order $1\lesssim Q\lesssim 10$ . Therefore, by analogy with the result for the uniform geometry, for geometrical wavelengths significantly smaller than these $O(1)$ thermal entry lengths, the localisation into features on the length scale of the topography is suppressed.

3.3. Long-wavelength regime

In the regime of long wavelength, $\lambda \gg 1$ , the different regions of the fissure are unable to ‘communicate’ with one another and a leading order analysis can be undertaken in which $x$ -derivatives vanish, and every vertical slice of the fissure can be considered as a planar problem of a different thickness. This regime is an interesting end member and could represent the variation of fissure width on the length scale set by the full length of the fissure – typically fissures are widest at their centres and get increasingly thin towards their ends. When the wavelength of the perturbation to the width is long, we can rescale the $x$ -coordinate by $\lambda$ and consider the variables to be asymptotic series in $1/\lambda \ll 1$ . Taking only the leading order of the governing equations, we obtain

(3.1) \begin{equation} \frac {\partial {\mathcal{E}}}{\partial t}+\Pi \frac {\partial }{\partial z}\left ( \mathcal{S}({\mathcal{E}})\frac {\partial p}{\partial z}\right )=\frac {2\Theta }{\delta }, \quad \frac {\partial }{\partial z}\left ({\mathcal{L}}({\mathcal{E}})\frac {\partial p}{\partial z}\right )=0. \end{equation}

Namely, the equations for the planar problem, with $x$ entering only as a parameter via the width, $h(x)$ . These equations and the resulting steady-state flux, $Q$ , can then be rescaled to the case with $h=1$ , by writing

(3.2) \begin{equation} {\mathcal{E}}=h\tilde {{\mathcal{E}}}, \,\,\delta =h\tilde {\delta }, \,\,{\mathcal{S}}=h^3\tilde {\mathcal{S}}, \,\,{\mathcal{L}}=h^3\tilde {\mathcal{L}}, \,\,\Pi =\tilde {\Pi }/h^4, \,t=h^2\tilde {t}, \text{ and } Q=\tilde {Q}/h. \end{equation}

Thus, given the uniform steady-state solution curve $(\Pi _0, Q_0)$ , we can find a steady-state solution curve for the position $x$ in the non-uniform fissure, as $(\Pi _0/h^4, Q_0/h)$ . In particular, at the widest region, we have $h=1+A$ and at the thinnest region, we have $h=1-A$ . Figure 7(a) shows the result of this rescaling for $A=0.02$ : the wider geometry (dashed) exhibits a higher flux and consequently transitions from the fast branch to the slow branch via a pair of fold bifurcations at a lower pressure drop; conversely, the narrower geometry (dotted) exhibits a lower flux and transitions between the branches at a higher pressure drop. If we then consider the average of the fluxes at each pressure drop, we obtain a new curve shown by the solid curve. The curve in figure 7(b) shows the prediction for the ratio of fluxes, $Q_r\equiv Q_+/Q_-$ , in the wider and thinner regions of the fissure, from this long-wavelength analysis. From these two panels, we observe that the solution curve exhibits the key features of the non-uniform geometries at larger wavelengths, namely a stable middle branch with significant flux localisation (large $Q_r$ ), connected to the fast and slow branches by two pairs of fold bifurcations. Also shown as points in each panel are the results of numerical solutions for $A=0.02$ and $\lambda =20$ . These were obtained by natural continuation in decreasing $\Pi$ , so the whole solution curve is not obtained; however, the branches that are obtained match the predictions well, closely capturing the maximum and minimum fluxes, $Q_+$ and $Q_-$ , and their ratio, $Q_r$ , and showing roughly the same critical pressure drops where the solutions change from one branch to another.

Figure 7.

(a) Steady-state solution curve (as in figure 3) for the uniform geometry, $h=1$ , re-scaled to the width of the widest, $h=1+A$ (black dashed), and thinnest, $h=1-A$ (black dotted), regions of the non-uniform geometry with $A=0.02$ and $\lambda \to \infty$ according to (3.2). These represent the maximum, $Q_+$ , and minimum, $Q_-$ , surface fluxes in the non-uniform geometry. The black solid curve shows the average of the fluxes, $(Q_++Q_-)/2$ . Points are shown from numerical solutions with $A=0.02$ and $\lambda =20$ , obtained by natural parameter continuation in decreasing $\Pi$ ( $Q_+$ , squares; $Q_-$ , crosses; average, circles). (b) Ratio of maximum and minimum flux, $Q_r=Q_+/Q_-$ , as a function of $\Pi$ . The line is obtained from taking the ratio of the fluxes represented by the dashed and dotted lines in panel (a), and the circles are obtained from the numerical simulations, as in panel (a). (c) Average flux per unit length, $Q$ , as a function of $\Pi$ . The lines are obtained by the long wavelength approximation (3.4). The different curves correspond to continuation in $\Pi$ in either the decreasing (solid) or increasing (dashed) directions. Also shown as circles are the numerical results as in panels (a) and (b). (d) Colour plot of scaled surface flux ratio, $\tilde {Q}_r$ (defined by (2.29)), as a function of $\Pi$ and $A$ for numerical solutions with $\lambda =20$ . The solid red lines show the boundaries predicted by the long-wavelength analysis ( $\Pi =241/(1+A)^4$ and $\Pi =241/(1-A)^4$ ), and the red dotted line shows the predicted lower boundary offset by a pressure drop of $5$ . Note the logarithmic colour scale. See table 2 for other parameters used in these solutions.

We note that the use of the average of the largest and smallest fluxes, $(Q_++Q_-)/2$ , in figure 7(a) is not the same as the average flux per unit length, $Q$ , used previously, but was chosen here for illustrative purposes. To obtain the prediction for $Q$ , we must account for the continuous variation of the width, $h$ . Suppose we have the flux as a function of the pressure drop, $Q_0(\Pi _0)$ , from the uniform geometry results (the choices made in defining such a function are discussed later). To obtain a long wavelength prediction for the average flux per unit length of the non-uniform geometry, we can then calculate

(3.3) \begin{equation} Q(\Pi ) \approx \frac {1}{L_x}\int _0^{L_x} \frac {Q_0(h(x)^4\Pi )}{h(x)}\, \textrm {d}x, \end{equation}

which, for the sinusoidal width, $h=1+A\cos (2\pi x/\lambda )$ , is independent of wavelength

(3.4) \begin{equation} Q(\Pi ) \approx \frac {2}{\lambda }\int _0^{\lambda /2} \frac {Q_0(h^4\Pi )}{h}\, \textrm {d}x=\frac {1}{\pi }\int _{1-A}^{1+A} \frac {Q_0(h^4\Pi )}{h\sqrt {A^2-(1-h)^2}} \,\textrm {d}h . \end{equation}

We note that the solution curve $(\Pi _0,Q_0)$ does not immediately provide a well-defined function $Q_0(\Pi _0)$ , due to the multiplicity of steady states. There is, however, a natural choice by which to define two such functions, corresponding to natural parameter continuation in either decreasing or increasing $\Pi$ . For the former, wherever $Q_0(\Pi _0)$ is multivalued, we choose to take $Q_0$ from the fast, hot branch, while in the latter, we instead evaluate the flux from the slow, cold branch. The resulting predictions for the average flux, $Q$ , for the sinusoidal geometry with $A=0.02$ are shown in figure 7(c). Here, the solid line is for the natural continuation in decreasing $\Pi$ and is compared against the fluxes obtained from the numerical solutions with $\lambda =20$ (circles, also obtained by continuation in decreasing $\Pi$ ). These show good agreement. The dashed line, for comparison, shows the prediction for natural continuation in increasing $\Pi$ , giving a lower flux where the two differ.

For the current parameters and the uniform geometry, $h=1$ , the fold bifurcation on the fast branch occurs at $\Pi \approx 241$ at which the flux drops from $Q\approx 6.6$ to $Q\approx 1$ (when continuing naturally in decreasing $\Pi$ ). Thus, for the long-wavelength approximation, we anticipate the first fold bifurcation to occur at $\Pi =241/(1-A)^4$ . At this bifurcation, the thinner region exhibits a reduction in flux by a factor of $6.6/(1-A)$ , while the wider region maintains the same flux. Thus, the surface flux ratio increases by a factor of $6.6/(1-A)$ when the solution drops onto the focussed branch at this point. The second fold bifurcation is predicted to occur at $\Pi =241/(1+A)^4$ , where the surface flux ratio decreases by a factor of $6.6/(1+A)$ . Thus, as the amplitude, $A$ , decreases towards zero, the stable focussed branch shortens but persists, ultimately shrinking to the point $\Pi =241$ , with $Q_+=6.6$ and $Q_-=1$ , as $A\to 0$ . The prediction for the critical values of $\Pi$ which bound the focussed branch is compared with numerical solutions for $\lambda =20$ in figure 7(c) showing good agreement for the upper boundary, $\Pi =241/(1-A)^4$ , but weaker agreement for the lower boundary where the long-wavelength analysis under-predicts the critical pressure drops observed in the numerical solutions by roughly $5$ at all amplitudes (this is also noticeable in the other panels of figure 7 for the specific choice of $A=0.02$ ). It is unclear what causes this discrepancy. However, since the thinner region lies on the slow branch of the solution curve, it becomes unstable to the thermoviscous fingering instability, resulting in finger structures within the cool, slow region of the fissure. The wavelength of these fingers decreases with decreasing pressure drop, $\Pi$ , and so it is possible that the assumption of large length scale in the $x$ -direction becomes invalid at the boundary between the hot and cold regions as the critical pressure drop is approached. In general, if gradients are greater in the $x$ -direction than implied by the rough scaling, $1/\lambda \ll 1$ , this will promote heat fluxes between the regions, resulting in a transition to the uniformly slow branch at larger $\Pi$ . Nonetheless, the predictions of the long-wavelength analysis give a decent agreement to the numerical solutions, particularly for the critical pressure drop resulting in the onset of localisation, which is perhaps of greatest concern for the volcanological application.

Figure 8.

(a,b) Scaled surface flux ratios, $\tilde {Q}_r$ , as a function of the dimensionless pressure drop, $\Pi$ , for steady and periodic solutions with $A=0.02$ and $\lambda =2$ over two ranges of $\Pi$ corresponding to the absence of stable steady states in figure 5(b). For periodic solutions, the maximum and minimum values of $\tilde {Q}_r$ are plotted as solid lines, and the average over the period is plotted as a dotted line. Markers indicate the bifurcations at which the steady branches lose stability, with triangles corresponding to infinite period bifurcations and the dot indicating a supercritical Hopf bifurcation. (c,d) Period of the solutions as a function of $\Pi$ , with the two panels corresponding to the branches shown in panel (a,b). Red dashed lines indicate the ends of the stable steady-state solution branches at either side. See table 2 for other parameters used in these solutions.

3.4. Periodic solutions

Here, we detail the origin and the bifurcation structure of the periodic solutions evidenced by the absence of steady states at points marked in red in figure 5(b). We take two examples of branches of periodic solutions along the solution curve for $A=0.02$ and $\lambda =2$ , tracking the solutions as $\Pi$ is varied. For these solutions, we took $L_x=1$ , i.e. the smallest domain that captures the $\lambda =2$ wavelength, given the symmetry boundary conditions at either end. The scaled surface flux ratio, $\tilde {Q}_r$ , and the period are plotted as functions of $\Pi$ in figure 8, showing that one periodic solution branch persists over the range $202.6\lesssim \Pi \lesssim 216$ and two different periodic branches occur over the range $229.4\lesssim \Pi \lesssim 233.8$ . At $\Pi \approx 233.8$ , on the steady solution branch, a pair of complex eigenvalues of the stability problem cross the imaginary axis, indicating a Hopf bifurcation. The amplitude of the stable periodic solutions obtained near $\Pi =233.8$ grows from zero, and their period starts at that of the unstable periodic mode at the bifurcation, $2\pi/\textrm{Im} (\sigma )\approx 38$ . Thus, the bifurcation to the periodic solution is a supercritical Hopf bifurcation (e.g. see Strogatz Reference Strogatz2015b ). The periodic solution branch arising at the Hopf bifurcation is short-lived, becoming unstable to a fold bifurcation at $\Pi \approx 233.45$ , resulting in the system transitioning to a new stable periodic solution branch, with a longer period and larger amplitude. At the other end of this periodic solution branch, the periodic solution is destroyed in an infinite period bifurcation at $\Pi \approx 229.4$ , where a saddle-node point appears on the periodic orbit and the system follows the new stable steady-state solution branch at lower $\Pi$ . The period-averaged scaled surface flux ratio, $\tilde {Q}_r$ , is naturally continuous at this bifurcation, since the periodic solution spends an increasing proportion of its (diverging) period close to the state that becomes the stable steady solution after the bifurcation. For the second periodic solution branch, shown in figure 8(a,c), the bifurcation structure is slightly different, with the transition between the steady-state and periodic solution branch resulting from an infinite period bifurcation at both ends.

Figure 9 shows the dimensionless average temperature, ${\mathcal{T}}(x,z)$ (see (2.13)), for a selection of solutions along the same solution branch (i.e. for the sinusoidal geometry with $A=0.02$ and $\lambda =2$ ). Panels (a)–(c) show steady solutions at $\Pi =200$ , $220$ and $234$ , representing solutions on the three steady-state branches separated by the periodic solution branches in figure 8. These indicate how the spatial structure of the fingered states changes after passing through the periodic solution branches. The other panels show three states obtained at particular times (with $t=0$ chosen arbitrarily) during the periodic solutions at $\Pi =210$ (panels d–f) and $\Pi =232$ (panels g–i), indicating that these solutions exhibit oscillations between states that are similar in structure to stable steady solutions at larger and smaller values of $\Pi$ . While the exact details of two fingered solutions might vary subtly, there are two main factors that govern the broad configuration of a solution. The first is the number of fingers in the domain or, equivalently, the average wavelength of the fingers. The second is the location of these fingers, which (given the symmetry boundary conditions) largely corresponds to a binary choice of whether $x=0$ corresponds to a (relatively) hot or cold region of the fissure. As for the uniform case (see figure 3 b), the preferred wavelength of fingers scales with the thermal entry length and, thus, the average wavelength increases as $\Pi$ increases. For example, at $\Pi =200$ (figure 9 a), the solution exhibits five hot fingers per wavelength of the fissure geometry (half of this wavelength, $\lambda =2$ , is shown in each panel), while at $\Pi =220$ and $\Pi =234$ (figure 9 b,c), there are four. In the non-uniform geometry, the preferred wavelength also depends on the local width of the channel, with longer wavelengths preferred in wider regions (since the thermal entry length also increases with channel width). This effect is visible in figure 9, where the wavelength of the fingering pattern generally gets shorter between the widest region of the channel, at $x=0$ , and the narrowest, at $x=1$ . The dependence on $\Pi$ of the choice between a hot or cold region at $x=0$ is less clear. Naively, one might anticipate that the most stable configuration would be one in which the widest point of the channel is hot and/or the narrowest point is cold. However, given the continuous variation of fissure width, the situation turns out to be more subtle than this. For example, at $\Pi =200$ , both choices result in stable steady states. Similarly, both configurations arise during the periodic solution at $\Pi =210$ (see figure 9 e,f). In figure 9(a), we show the state with a cold region at $x=0$ , which is the solution reached when continuing in $\Pi$ from the periodic solution branch at $202.6\lessapprox \Pi \lessapprox 216$ (i.e. it is a state similar to that shown in figure 9(f), and not panel (e), that becomes the saddle-node point in the infinite period bifurcation at $\Pi \approx 202.6$ ). This is also the configuration which is the more stable of the two at $\Pi =200$ (in the sense that the growth rate of its least stable mode is more negative). In general, the preferred configuration depends on $\Pi$ . For example, both the solutions at $\Pi =220$ and $\Pi =234$ in figure 9(b,c) have an average wavelength of 0.5, but the solutions have opposite configurations. There are thus several (potentially) competing influences on the configuration of fingers in a flow-focussed state: the preferred wavelength of fingers in an average sense; the preferred wavelength of fingers locally, due to the local channel width; and the relative stability of the two choices of finger location. The periodic solution branches then arise where these influences are in conflict and no choice of configuration is stable, with the solution instead oscillating between multiple, weakly unstable states. During a periodic orbit, the evolution between these states happens in two main ways: a propagation of the fingers down the length of the fissure (e.g. panels e ,f and g, i in figure 9); or the removal or addition of a finger at one or both ends of the domain (e.g. the evolution to and from panels d and h). As demonstrated in figure 8, the periodic solution branches can evidence a range of different behaviours, and we do not claim to have provided an exhaustive exploration of these states. We also note that no periodic solutions were observed for the uniform geometry. Two of the factors previously discussed (the spatially varying preferred wavelength, and the choice for the location of the hot and cold fingers) are absent without the presence of wider and narrower regions of the fissure, and so it is perhaps unsurprising that the uniform geometry does not exhibit oscillatory dynamics.

Figure 9.

Colour plots and contours of average temperature, ${\mathcal{T}}$ , for several solutions in the sinusoidal geometry with $A=0.02$ and $\lambda =2$ . (a–c) Steady solutions at (a) $\Pi =200$ , (b) $\Pi =220$ and (c) $\Pi =234$ . Panels (d)–(f) and (g)–(i) show snapshots at several times in the periodic solutions at $\Pi =210$ (period $\approx 182$ ) and $\Pi =232$ (period $\approx 115$ ), as indicated in the titles, and $t=0$ chosen without loss of generality. Contours are plotted from ${\mathcal{T}}=0.1$ nearest the top of each panel, increasing in increments of ${\mathcal{T}}=0.1$ . The axes are the same in all panels, and are omitted from all but panel (a) to avoid crowding. See table 2 for other parameters used in these solutions.