Fluid Flow With Three Upstream Configurations in Freezing Tubes

The accumulation of frozen liquid around a central passageway of melt as it flows through a freezing region can make calculations very challenging. To both illustrate and to quantify some of these challenges from freezing, a model equation is developed. It simplifies the solution of Holmes (2007, https://gfd.whoi.edu/wp-content/uploads/sites/18/2018/03/MHolmesGFDReport_30151.pdf) for low Reynolds number single component liquid flow through a long tube that has a wall kept at subfreezing temperature. This model equation is used in conjunction with three different upstream configurations, each with parameters expressing their behavior. Analytical and numerical results give the parameters that have criteria for: the freezing of a compressible upstream reservoir that includes oscillatory behavior; the freezing of flow fed through a constriction with a large upstream pressure, just like a dripping water faucet during winter; the evolution of flow in multiple tubes connected by an upstream manifold, where some tubes end up with full flow and others freeze shut. Numerical runs with 1,000 tubes give a formula for the spacing between actively flowing (non‐frozen) tubes over wide ranges of the two upstream parameters (flow rate and manifold resistance). Results have implications in various areas in earth science. Some are: oscillatory and freezing shut criteria for flow of magma from a compressible region, a criterion for wintertime ice accumulation at natural springs, and the spacing between volcanos.

Obviously, freezing must constrict the flow region and increase resistance. The same increase in resistance occurs if, instead of freezing, the fluid has a viscosity that increases with colder temperature. In that case, for progressive cooling, a flow becomes focused into narrow channels surrounded by colder, more viscous sluggish flow. Therefore, the focus into more constricted regions is similar to the focusing by solidification. Various geometries that have been studied of fluid flow with temperature-dependent viscosity include: regular circular slots (Helfrich, 1995;Whitehead & Helfrich, 1991;Wylie, Helfrich, et al., 1999;Wylie & Lister, 1995); gelatin (Pansino et al., 2019 and citations therein); and cracks . Therefore, flow with viscosity variation with temperature can also be sensitive to upstream dynamics.

A Laboratory Demonstration of Channeling
A laboratory experiment for teaching exhibits transition from wide liquid flow to a flowing channel surrounded by solid. A wax is injected with a positive displacement pump at a constant rate into the center of a circular slot over a carefully leveled aluminum disk 0.4 m in diameter painted black and kept at a temperature below the solidus (Figure 1a). The slot (of fixed small thickness approximately 2 mm) is between a transparent circular polycarbonate lid and the disk. The layout is similar to previous experiments with paraffin (Whitehead & Helfrich, 1991) and flow of oversaturated water (Kelemen et al., 1995). Both demonstrate the formation of a channel. The liquid is forced to spread from the center outward over the cold disk and to make its way to the outer edge. Most of the liquid solidifies, but at least some of it flows all the way to the edge and spills into a catch basin. The laboratory liquid, 1-hexadecene, is a clear liquid at room temperature and becomes a white waxy solid at 3.6°C. The cold disk is at −5°C. The volume flux rate is 9.1 × 10 −6 m 3 s −1 . Figure 1 shows the sequence of liquid flow and solid accumulation. After the pump is started, one frozen fan of wax accumulates (Figure 1b), followed by a new outbreak of flowing liquid leading to a second fan ( Figure 1c). Then, there are many subsequent cycles of outbreak-fan formation (Figure 1d) so that the sequence of fan formation and outbreak ultimately circles around 360°. After 45 min, the total region ends up being filled with solid. At that point, flowing liquid occupies a comma-shaped region near the center that was laid down during the fan sequence. The video (in Supporting Information) shows that next, the liquid forces the lid upward a small amount because the positive displacement pump feeding in the melt at the center can produce immense pressure when all the material is frozen. A very thin gap between the solid and lid opens and an axisymmetric flow of melt goes radially outward (Figure 1e). This radial flow is almost immediately followed by the appearance of one rapidly amplified dark drainage channel extending from the central hole to the outside rim of the cylinder. The channel becomes progressively darker and wider during WHITEHEAD 10.1029/2020JF005969 2 of 25 a five-minute period as the channel melts its way down through the wax all the way to the aluminum disk. Thereafter, the entire flow occupies this channel of fixed width (Figure 1f). Additional runs have a width of the final channel proportional to flux rate (C. J. Mills, private communication). Similar results are described in Kelemen et al. (1995) with ammonium chloride.

This Study
The purpose of this study is to illustrate the influence on freezing flow by three different upstream regions. Generally, the flow in freezing regions is difficult to calculate, so to simplify the freezing dynamics, the simplest geometry, tube flow is used along with simplified mathematics. The freezing flow solution in Holmes (2007) in a tube is replaced in Section 2 by a model of flow that replaces the complicated calculation of the thermal fields with analytic functions. The Appendix shows the analysis by Holmes (2007) that leads to our simple model. Then, three sections show analysis of freezing flow with different upstream configurations. Section 3 analyzes the stability properties of this model when it is fed by an upstream storage chamber with a free surface in the field of gravity. This is equivalent to a compressible reservoir like the upstream WHITEHEAD 10.1029/2020JF005969 3 of 25  Holmes-Cerfon & Whitehead (2011). When performing numerical calculations with the model equations, freezing shut leads to pressure approaching infinity. There are situations where this causes difficulty that is overcome by either terminating the calculation and setting flow to zero with complete freezing, or by continuing the calculation and adding another physical process into the model to avoid the large pressure. We added the physical process of not allowing the radius for the flowing liquid to be smaller than a "minimum radius." The flow with this radius is seepage flow that has very tiny volume flux compared to the other flows in question. This helps mathematically because as flux rate becomes very small, the minimum radius causes the pressure-flux rate curve to bend down and approach zero instead of extending up to infinity. Thereby, seepage flow allows calculations to continue forever. Numerical calculations produce oscillations like those with viscosity-temperature variation in the laboratory (Whitehead & Helfrich, 1991). Section 4 has a second upstream configuration that is like a dripping faucet in freezing weather. The criterion for freezing up/seepage flow is found and explained. Section 5 analyzes flow and freezing up/seepage flow for multiple tubes (from 2 up to 10 4 ). These are aligned next to each other and fed by a manifold that connects them together in the upstream region. In these calculations, the minimum radius and consequent seepage flow avoids pressure in the manifold going to infinity when flow is approaching total freezing, which produces cross-manifold flow rates going to infinity. The numerical results produce a formula relating the spacing of active tubes to the parameter expressing a resistance coefficient of each manifold tube divided by the upstream volume flux rate. Results are applied to some problems in igneous flow.

Previous Solutions
The model is a simplification of one of the simplest examples: a liquid flowing through a pipe held below the liquid solidus temperature. The mathematical solutions for this configuration were first developed by Zerkle and Sunderland (1968), Sakimoto and Zuber (1998), and references therein based on separation of variables with eigenvalues and eigenfunctions by Graetz (1883).
The analysis is valid for Peclet number   0 / Pe ur of order one and Prandtl number Pr >>1, consequently,

Re ur
Here, u is velocity,  is kinematic viscosity,  is thermal diffusivity, and 0 r is tube radius.
The developing flow at the tube entrance is also ignored so Graetz number    2 0 / 1 Gr ur L with L the tube length. The same limits apply throughout this paper and in addition all fluid properties are constant and independent of temperature. Holmes (2007) and subsequently Holmes-Cerfon and Whitehead (2011) calculated the flow in these limits into a freezing tube with constant viscosity. Accumulation of solid produces a decrease in fluid-solid radius in the flow direction ( Figure 2). The central attribute that leads to instability of these flows is a pressure minimum at some value of flux rate with pressure p approaching infinity as volume flux rate approaches the two limits of  0, as in Figure 2c; large flow rate makes large pressure (pipe flow) and tiny flow rate makes tiny liquid radius and large pressure. The variables in Figures 2b and 2c T the temperature of the tube wall, and T i the temperature of the inflowing liquid.

The Simplified Model
For this study, we adopt a simplified formula with one curve replacing the curves like those for different T n in Figure 2c. Instead of an annulus of solid that has a decreasing inner radius in the flow direction, the radius varies only in time and not along the tube. The dimensionless equivalent of Equation A2 has a balance WHITEHEAD 10.1029/2020JF005969 4 of 25 between volume flux q, the radius of the liquid-melt interface a and the pressure drop across the tube p is independent of the flow direction. It is The dimensionless Stefan condition (see Equation A6) for the evolution of the fluid/solid interface is a balance between the growth or decay of a solid with latent heat of solidification and the divergence of the conductive heat flow there (see also Turcotte & Schubert, 2002, p. 162) The time is scaled using L H is latent heat of solidification and C p is specific heat. Equal values of specific heat for solid and liquid are used for simplicity. Also, E is dimensionless radial heat flow in the solid and I is dimensionless radial heat flow in the liquid at radius a. With Stefan Number small, the radius a, which is at the melting temperature, changes at a time scale slower than the thermal conduction timescale so that a steady heat flow occurs in the solid along the tube as in Holmes (2007) and Holmes-Cerfon and Whitehead (2011). In the solid, in the Appendix) is replaced by a function that does not vary along the flow direction. The first term in a Taylor series expansion about ln(a) produces The relation is best for a close to 1 with the values changing significantly from Equation A8 as  0 a .2. Finally, the term   1 / a t in front of E is also set to 1. Physically, this means that the heat flow formula governs the heat flow in cartesian coordinates through a slab over the inside area of the tube with a liquid-solid radius close to r = 1. It also means that heat flow (cooling of the flowing liquid) is larger than would have occurred otherwise. Next, for the heat flow balance in the liquid, hot liquid flows in from the left and leaves at a lower temperature on the right. The greatest possible heat loss is with     , / 4 I q q in Equation 1 with the liquid exiting at the solidus temperature. Therefore, Equation 2 is simplified to WHITEHEAD 10.1029/2020JF005969 5 of 25 Using the rescaled variables   4 n q q T ,   4 n p p T , and   n t T t, Equation 1 is the same form but with primes but Equation 3 becomes (4) As emphasized above, the compelling justification for these simplifications is pragmatic rather than physical. The approximations are clearly not rigorous. Our best justification is that Equation 3 has steady flows that produce a realistic q-p curve that is easily used for stability studies. The curve using Equations 2 and 3 ( Figure 3c) is a small distance to the left of the curve with  0.1 n T .
The fundamental objective of this study is to use Equation 1 with primes and Equation 4 to explore the dynamics with the three different upstream configurations sketched in Figure 4. The first is a compressible storage reservoir lying upstream of the tube. The second is a fixed resistance in series with the tube fed by a reservoir at constant pressure. The third has multiple tubes connected by a manifold.

Compressible Upstream
The addition of a compressible upstream reservoir can be considered to be a model of a magma delivery system in the earth, and possibly to planets and moons, too. Time-dependence is a fundamental feature of magma production in the earth irrespective of composition, temperature and geometry. Many mechanisms can lead to time-variability in a system with steady forcing such as volatile content and outgassing, brittle behavior, viscosity variation, and crystal settling, but this model produces time dependence without them. Additional features such as outgassing and   Journal of Geophysical Research: Earth Surface viscosity variation might be added later to produce highly eruptive cycles with faster time scales (Wylie, Voight, & Whitehead, 1999).
The simplest upstream condition consists of a reservoir of fluid of dimensional area A with a free surface that can go up and down (it is essentially a compressible reservoir even though the fluid is incompressible, see Figure 4a). It is fed by a constant inflow u q , which is divided by 4T n (The prime in the definition of u q is left out to be consistent with a steady flow notation of the basic flow that is introduced in Section 3.1.) Fluid flows out of the reservoir and into the tube with volume flux rate  q . The pressure in the upstream reservoir obeys The dimensionless growth rate is

Stability with Compressible Upstream
Flow rate, radius, and pressure are expanded in a series expansion for their amplitudes with a steady component and a time-dependent component and from Equation 1 The shape of Equation 9 has the desired form shown in Figure 3c. The asymptotic log-log slope at large 0 q corresponds to simple tube flow independent of T n and is identical to all the solutions where the flow is so rapid that and Equation 5 is with roots (using Equations 6 and 9) The flow is unstable if growth rate  has a positive real part. Because  and 0 a are positive, the term is negative and the term under the radical sign has smaller real magnitude than the term to the left of the radical sign. If the term to the left is positive, then positive growth rate exists so the formula for instability is This indicates that if  is sufficiently large there is no instability. In addition, at zero growth rate, the radical term in Equation 14 is imaginary and the neutrally stable flow oscillates, although for sufficiently small  the growing instabilities are overdamped. Equation 15 can be rewritten using Equations 6-9 as Curves for five values of  c for ranges of u q are plotted in Figure 5 and the intersection of the curves with Equation 9 determine the minimum value of u q for stable flow for each value of  . c Positive growth, which leads to instability and presumably ultimately freezing shut, lies to the left of the intersection and stability to the right.

Journal of Geophysical Research: Earth Surface
The curve in Figure 5 is similar to the one in Figure 6a of Holmes-Cerfon and Whitehead (2011), who used the same upstream condition. First, flow is only possible above a specific inflow, here it is for q u > 3. Second, the dashed curves are very similar. Third, the limits are the same too. The limit   0 has pressure constant for all time according to Equation 5, and the flow is stable only with q 0 > 3 and unstable otherwise. On the other hand, the limit    has all flow rates stable. Flow rates equal q u according to Equation 5 and solutions occupy the entire curve. Therefore, the model equation has flow and stability ranges that are similar, but that disagree somewhat quantitatively with the results for the complete solutions in Holmes-Cerfon and Whitehead (2011). This helps to justify the use of this simple model equation instead of the full solution.

Numerical Results
Equations 1 and 5 are easily integrated ahead in time with forward finite-differencing and Equation 4 is used to calculate a new value of  q . The code validation results are given in the Supporting Information. One example, typical of many is shown in Figure 6. The oscillation amplitude initially increases. At t = 191.226, when amplitude becomes sufficiently large, there is an abrupt decrease in radius signifying a collapse toward freezing shut. This starts at the instant when the smallest radius occurs in the cycle. Although this figure is typical, some cases can be highly damped with perturbations decaying exponentially from the beginning. In all cases the sudden decrease signifies freezing shut and radius plunges toward zero.
At the freezing stage, the numerical calculation has radius shrinking to zero at a finite time ( Figure 6b). This is clear from Equation 4 because pressure remains finite from Equation 5 and therefore using Equation 1 the formula    3 / q a p a becomes zero and     / 1 a t . Therefore, at one particular time step the radius jumps to a negative value. In every one of our early numerical calculations, the calculation failed because Equation 4 crossed zero. This occurred no matter how small the time step was or no matter what the numerical method was for stepping ahead in time. At the freezing time step, there are two options. One is to set the flow value in the that tube to zero. In this case, before that happens, pressure continued to build up and the smaller the time step, the greater the pressure increase. Since a more accurate numerical code experiences a greater spike in pressure before freezing, this option is clearly unacceptable. The second option is to introduce additional dynamics. I decided to substitute a steady small minimum radius at every time step where the value of the radius becomes either negative or smaller than a fixed value. Substituting such a minimum radius always produces a small seepage flow that generates interesting new behavior without numerical failure. For the example in Figure 6, the minimum radius was first invoked at t = 191.226. After this, seepage flow continues and Equation 5 leads to a gradual increase in pressure (Figures 6c and 6d) that occurs until flow rate is great enough for the seepage flow to melt back and open the tube following Equation 4. This in turn causes the flow to become periodic because the minimum radius adds an additional straight line in the pressure-flux rate curve that extends from zero up to a point where it intersects Equation 9 as shown in the inset within Figure 5. After intersection, a new limit cycle oscillation occurs (Figures 6d-6f) with pulses of rapid flow separated by very slow flow. Figure 6f shows that the upstream pressure during the limit cycle is much greater than the original pressure. Numerous additional calculations showed us that the new limit cycle occurs throughout a wide range of parameter space.
The period of the limit cycle is not the same as the period of the linear instability. Instead, it depends on the minimum radius value so that the minimum radius is an additional parameter of the model. All aspects of this limit cycle are affected by minimum radius value including the time for build-up to the start of the limit cycle, the value of upstream pressure that is needed before the limit cycle begins, the limit cycle WHITEHEAD 10.1029/2020JF005969 9 of 25

Journal of Geophysical Research: Earth Surface
frequency, and the minimum and maximum values of flow rate and pressure for the limit cycle. The limit cycle involves a melt-back of the solid when pressure builds up enough to make the seepage flow rapid enough. Surprisingly, this flow rate is less than the flow rate at the instant of the beginning of the freezing shut event. This is clearly shown in Figure 6e. Apparently the flow rate at melt back occurs when the linear flux versus pressure curve for the minimum radius intersects the far left end of the curve for steady flow as sketched in the inset in Figure 5. This aspect is noted also by Helfrich (1995) for flow focusing with temperature-dependent viscosity.
The cycles are similar to oscillations in tube flow with temperature-dependent viscosity and upstream compressibility (Figure 7). There, instead of a minimum radius and seepage flow, the slow flow produces a cold very viscous "plug." This flow has a smooth p-q curve without discontinuous slopes that the cusp from the intersection of a straight line and Equation 9 has our model, sketched in the inset of Figure 5. Both of them seem to produce the same behavior. A similar increase in flow resistance occurs both with supercooled ice that produces crystals (Gilpin, 1981) and in lava flow situations with crystal accumulation and volatile excretion (Wylie, Voight, & Whitehead, 1999).

Formulation
The second upstream condition imposed here has the configuration in Figure 4b. It is inspired by the very well-known flow of water in pipes and in natural springs that persists during freezing temperatures. In fact, a common trick used by homeowners and plumbers to prevent pipe rupture during periods of freezing is to leave a water faucet with a dripping rate that is quite small for small ranges of subfreezing temperature or for short durations so that the water in the pipe does not freeze shut. In another example of a similar process, water continues to flow out of rock fractures long after air temperatures fall to below freezing, resulting in large accumulations of ice. These can become hazards in subfreezing railroad and highway road cuts, with some of them reaching great size. A hint of why flow exists with below freezing temperature is found in the limit of large  (Section 3) which is equivalent to an imposed steady flux rate where flow continues for any value (Epstein & Chueng, 1983;Holmes, 2007;Holmes-Cerfon & Whitehead, 2011). Therefore, an analysis of this problem that includes upstream dynamics of the dripping water pipe is useful.
Second, the upstream constriction, representing the valve in a faucet, can be pictured as a tube of radius f r and length f L . The dimensionless faucet pressure drop is thus The freezing tube and the faucet (either upstream or downstream) are connected in series to a reservoir at fixed large upstream pressure  u p so that

Stability
In general, the straight line has two intersections over a wide range of p u , one intersection at a tangential point and no intersections over the rest of the range of p u . For stability, Equations 10 and 11 for the small perturbations are used along with Setting   1 1 , t q p e , and combining Equations 10, 11 and 21 and then and using Equations 6-25 to simplify the coefficients, the formula for growth rate is Because the slope of Equation 9 is growth rate is simplified to This equation has a simple physical interpretation. Simply start at the zero flux axis with the line    Instability does not necessarily imply freezing. The perturbation 1 p has opposite sign from 1 q so that a perturbation with greater flux rate will have pressure decrease and tend to follow Equation 9 around to meltback and the steady flow. However, a perturbation with negative flux rate will have pressure increase and follow Equation 9 to smaller flows and thereby tend to freeze shut. Numerical calculations verify the values in Figure 8b.
Summarizing this section, faucet resistance R is quantified by a virtual radius with freezing shut very sensitive to its value. At sufficiently large values of upstream pressure p u for fixed R, the freezing is prevented by the same mechanism as with upstream pressure conditions alone. As p u decreases, the portion of pressure drop through the faucet (Figure 8a) begins to decrease and the pressure drop p 0 through the tube increases because a decreases. Finally, p u reaches a value that is small enough for the pressure drop in the tube to reach to its extremum for that value of R. Up to this point, the flow has been stable and small perturbations decay in time, but here Equation 24 indicates that the perturbation reaches neutral stability. Below this upstream pressure, no steady flow is possible and in addition the perturbation grows. This leads to smaller and smaller flow until presumably freezing occurs. Therefore, every faucet setting, given by a value of R has a critical value of p u .
For a faucet, its equivalent radius compared to tube radius could be much less than 10 −2 , resulting in R > O(10 8 ). Hence, flow freezes shut when the pressure and flow rate are reduced to the point with a large negative slope on the left-hand branch, which Figure 8 shows is considerably to the left of the minimum. Freezing shut also occurs if the initial steady flow is small enough to lie to the left of R c .

Multiple Tubes
Branching tubes are a model of sheet flow. Holmes (2007), numerically calculated flow in branching tubes where the source was comprised of a manifold connecting a large number of tubes. The manifold was simply tubes at the upstream temperature connecting the upstream tubes together. It received uniform inflow along its entire length. The mathematical solutions were numerically stepped ahead in time to see the evolution of flows. Fifty identical tubes responded subjected to influx values that were small enough so that only six or seven active tubes with volume flux rates to the right of the minimum in Equation 1 resulted with the rest freezing up. The calculations verified the expectation. It was necessary to set to zero the flux of any tubes that were freezing up and letting the pressure distribution along the manifold be determined by active tubes alone. Helfrich (1995) calculated planer flow with fluid having viscosity variation. This achieved flow focusing into discrete locations. Both results suggested that flow focusing is a topic that could be fruitfully quantified over wider ranges of parameters since only a small number of cases were considered. They motivated this study of multiple tubes connected by a manifold.

Two Tubes-Formulation and Analytical Results
Consider two tubes each fed by a source with flux rate q u with their upstream ends connected together by a "manifold tube" with flow back and forth (Figure 4c). An upstream pressure at the source, although possibly interesting, has not been analyzed yet. The relations corresponding to primed Equation 1 are pressure  1 p for tube 1, and, using similar notation for the pressure in tube 2, and also for the flux rates and radii of both tubes.
The manifold tube is kept at the upstream temperature and has different length and radius than the cooled tubes. Manifold flow resistance is inversely proportional to a resistance coefficient defined as C =  4 4 0 / m m L L r with  m the dimensional radius of the manifold tube and L m the physical length of the manifold tube. The two upstream conditions are Expanding as before, the notation of the previous sections can be used by adding a second subscript for the ith tube so for example    Examples for four values of C are shown in Figure 9. Intersections lie above the minimum and it is simply the inverse of the slope Equation 23. Some values of C have interesting behavior. First, the limit of large C is a horizontal straight line with two steady solutions, obviously only valid above the minimum so that in this limit

Journal of Geophysical Research: Earth Surface
the minimum radius has not been discussed so far in this section, this "impossible" dilemma is resolved by adopting the minimum radius that results in one large flow in partnership with one seepage flow. Consequently, the minimum radius is used in all the numerical calculations of multiple tubes.
The range of possible steady flows has been found, but are they stable? The steady flows have equations governing small time dependent perturbations that are first, the equivalents of Equation 10 for each tube (i = 1,2) and second the equivalent to Equation 11 The conditions in the upstream tube connecting them are It is convenient to modify Equation 37 using the equivalent of (2.5) to eliminate 0i a  

Journal of Geophysical Research: Earth Surface
Unfortunately, the algebra for two unequal flows is very complicated and is not developed but the analysis for two equal flows is straightforward.
The growth for this radius difference is positive if the value within the square brackets is negative, which becomes, after some manipulation a familiar formula Rewriting this using Equation 31, positive growth for instability requires which has a margin at Equation 36, with selected values shown by dashed curves in Figure 9. Their intersection with the steady flow curve (bold) gives values of the critical flow rate and Equation 44 shows that this occurs exactly at the tangent to the curve. Therefore, for both two tubes with identical flow rates and for the dripping faucet, a steady flow is stable in the entire range where the upstream volume flux rate is large enough to satisfy the steady flow equations.

Two Tubes, Numerical Results
The numerical calculation advances the two values of a by one time step using Equations 27 and 28 using forward finite-difference and then calculates q using these formulas derived from Equations 25, 26, and 30.
Then, the new values determine both pressures at the new time. In practice, one tube might begin to freeze and end up with radius shrinking rapidly toward zero when seepage flow occurs with the minimum radius (See inset in Figure 9a) as in Section 3. The code is validated (see Supporting Material) by comparing the parameters for stability of flows with Equation 45. For q u = 2 the instability occurs with C > 16/27 = 0.59259. To start, the initial flows in the two tubes are specified and the value of critical C is numerically calculated by trial and error. Results show that stability depends on the initial flow values. Starting with   1 1 q and   2 3 q , instability occurs (up to the fifth decimal) for C > 0.30682; then with   1 1.8 q and   2 2.2 q , C > 0.57826; next, with   1 1.98 q and   2 2.02 q , C > 0.59245; and with  1.998, and 2.002, C > 0.59259. Therefore, the stability criterion depends on  1 q ,  2 q and C. A comparison is also made between stability prediction Equation 45 and numerical runs with values of minimum radius of 10 −3 , 10 −4 (the value used WHITEHEAD 10.1029/2020JF005969 15 of 25 in subsequent calculations) 10 −5 and even 10 −13 . They all agreed. The value of minimum radius does not determine stability.
Other numerical results over a wide number of parameters verify the analytic formulas in Section 5.1. The Matlab code is quite simple and written out in the Supporting material. Figures 9b-9d has examples for three parameter pairs: one starts with unequal flows that approach balanced flows (panel c); and two cases start with almost identical flux rates that go unstable to having all the flow in one tube and seepage flow in the other tube (panels b and d). Finally, in no case have two unequal flows that satisfy the straight line intersections in Figure 9a remained steady. Flows always evolve to either two equal flows, or to full flow in one tube and seepage flow in the other.

Many Tubes-Numerical Results
Numerical calculations are easily formulated for more than two tubes. Each tube radius is advanced in time based on the radius and flux rate within each tube using equivalents of Equation 27 This resets flux rate for each tube after which the cycle is repeated.
To begin a numerical calculation, a fixed value of q u and C is specified and the initial radius for tube number i is given the flux rate q u (0.9995 + 0.0001var(i)) where var(i) is a random integer between 0 and 10 produced with a numerical random number generator. Radii and flux rates in each tube thereafter advance in time until steady state is reached. Validation of the program (in Supporting Information) was made by closely following each time step with a 10 tube manifold and checking that each individual tube flow obeys Equations 25-30.
When instability develops with some tubes having larger flows and others smaller ones, Equation 50 proceeds without interruption even after seepage flow develops. Invoking a minimum radius is essential since otherwise Equation 50 develops "shrinking denominators" as some radii become very small with a consequential immense increase in pressure leading to unphysically large manifold flows and numerical instability. Figures 10a and 10b shows a typical evolution of flux rate and radius for 101 tubes. The time step is small enough to allow different wavelengths of a perturbation to test for different growth rates as found with temperature-dependent viscosity (Helfrich, 1995;Wylie & Lister, 1995), but they are the same. All calculations exhibited no selective wavelength. Instead, the random perturbation profiles of both the flux rates and radii remain almost perfectly preserved during instability growth. The preservation continues throughout an "early stage" (up to t = 0.3 in Figure 10). This stage terminates at different times depending on perturbation WHITEHEAD 10.1029/2020JF005969 16 of 25 size, q u and C. Suddenly, at a time that depends on the parameters (from t = 0.35 to 0.4 in Figure 10), there is an "intermediate stage" where the profiles and radii have order one variation and they begin to dramatically change. Some radii and flow rates plunge toward zero and others increase. The evolution of each individual tube is not understood. For example, a tube radius might first decrease and then increase or vice versa as the upstream manifold pressure distribution readjusts. Last, Figure 10 shows that a late stage follows and at t = 2 some tubes approach seepage and others fully flow. Ultimately, all flux rates and radii in the active tubes become almost exactly equal and all seepage flux rates do too. A cross-manifold flow remains that distributes fluid from the uniformly spaced sources to the active tubes. In this model, the ends of a manifold have zero lateral flux rate and this exerts some influence not yet documented or understood. In spite of this, results are clear. For example, the steady final distribution at t = 2 for q u = 0.1, C = 1 ends with in flow in six tubes (Figures 10a and 10b). In smaller q u , only one tube has flow, and this persists even for q u small enough for seepage flow in the final tube. A sequence with an unchanging distribution in the early stages of a numerical model with viscosity variation of cylindrical-slab flow seems to be similar to this (Figure 14 of Helfrich, 1995) ending with one flowing region and where everything else is decaying away.
Our small minimum radius value of 0.0001 that is used for all calculations makes reproducible results that are subject, of course, to the limits of random initial conditions. The seepage flux rates are very tiny in the volume flux budget at the end of all calculations. For example, even for an extreme case with 1,000 tubes where only one tube remains active at the end of the freezing up/seepage flow sequence, more than 99% goes through the active tube and less than 1% of the imposed flux going through the 999 seeping tubes.
This evolution of small perturbations that evolve to flow that is equal in selected tubes with seepage flow in the others occurs in all 1,183 numerical runs. The results in Figure 10 are listed in the Supporting Tables. All runs span wide ranges of q u (10 −5 -100), C (10 −4 -10 8 ), and N (2-10,000). Each realization follows the nonlinear evolution ending in a few actively flowing tubes with equal flux rates ( Figure 10a) and radii (Figure 10b) at t = 2. Some effects of the manifold ends exist. A a variation of the spacing away from the center becomes more than 10% for N  1,000. With random initial conditions, the final number of active tubes, #, has a statistical spread in 1,000 realizations. This is shown for one case in Figure 10c. The result is insufficient to determine whether it is bell shaped, which might not happen because of the nonlinear evolution.

Journal of Geophysical Research: Earth Surface
Although flux rates and pressure of active tubes arrive at one point on the    p q curve, no rule is known governing the final values. For example, the rate in each tube in Figure 10a is about   1.7 q , which is below the stability value for two tubes in Figure 9a and thereby lies on the unstable branch. Therefore, the concept of "some flows grow in the stable branch and others decay on the unstable branch" does not hold. Helfrich (1995) also reports this for flows with fingering due to temperature-dependent viscosity.
After some searching, a systematic empirical dependence between tube number # (consequently spacing N/#), and the parameter group q u /C 1/4 was found for wide ranges of C (10 8 ) and q u (1.5 × 10 4 ). All tubes flow for q u /C 1/4 > 0.55 and flow fills fewer tubes for the remaining 154 runs. The trends in log-log space are linear, parallel and logarithmically close to linearly proportional to the parameter group q u /C 1/4 (Figures 10b  and 10d) in spite of no averaging as in Figure 10c for randomness.
Since C is defined to be proportional to 4 m a , C 1/4 will be called the "scaled manifold radius." The linear trends in Figures 10d and 10f have slopes proportional to   1 1/4 / u q C and active tube spacing N/# is linearly proportional to scaled manifold radius. To quantify the results further it is useful to note that each volume flux rate is simply q = Nq u /#. All radii are also equal so that the radius a for steady flow in each active tube is readily calculated using Equation 7. The ratios of this radius compared to the scaled manifold radius 1/4 / a C for the points shown in Figure 10 are shown in Figure 11a. The ratios are not constant, but they all are clearly of order one. For  1 C the ratio 1/4 / a C has considerable variation of a little over 2 with a total range from 0.3 to 0.68. For C = 100, the mean ratio is 0.225 with a standard deviation of 0.003. For C = 10 4 , the mean ratio is 0.082 with standard deviation 0.0088. Therefore, to a first approximation the radius within a flowing tube is linearly proportional to the scaled manifold radius C 1/4 with a proportionality constant (Figure 11a) that is order one.
The number of flowing tubes is not only statistical, but it is also influenced by history. Figure 11b shows contours of radius for all active tubes with q u set to four progressively lower values and it is a good illustration of the evolution of active tubes. When this run is continued with the opposite sequential increases up in q u there is hysteresis with no increase in the number of active tubes. This hysteresis is explained by considering that for flow in a single tube, the total flux rate is 0.04 × 200 = 8 making an upstream pressure of about 12 (see Figure 5). This pressure makes only a tiny seepage flux rate of 1.2 × 10 −15 but the seepage flux rate needed for the straight line of seepage flow to intersect Equation 9 exceeds 1.

Journal of Geophysical Research: Earth Surface
In summary of this section, the model curve in Figure 5 causes problems by building up huge pressure in the manifold. This produces instability characterized by immense flows back and forth. The imposition of a minimum radius ( Figure 5 inset) removes this difficulty. Consequently, the numerical calculations work well at documenting evolution for  1, 000 N . For q u /C 1/4 < 0.55 both the spacing between active tubes and the value of the active tube radius depends primarily on the scaled manifold radius C 1/4 divided by q u . For growth from random noise, the relation between q u and # has statistical results that cluster around a central peak.

Discussion
A simple model is used to analyze a number of flows with three different upstream conditions. Explicit formulas for stability and other aspects of each flow lead to insight into the dynamics of freezing for each of them. For a compressible upstream chamber, the two limits of imposed upstream pressure for small  and constant flux rate for large  are recovered. There is a range where instabilities oscillate similar to Holmes-Cerfon and Whitehead (2011). For the frozen water faucet configuration, freezing shut occurs when the pressure change with volume flux rate of the flow equals slope of the curve as shown in Figure 8. For branching tubes, Equation 25 indicates that freezing shut of one of a pair of tubes occurs if the inverse of the resistance coefficient between the two tubes upstream is greater than the tangential slope as in Figure 9. For all three configurations, numerical calculation for this model with finite time steps does not extend all the way to perfect freezing unless a special numerical addition is implemented to remove high pressures for very small flow. A minimum radius is used to allow numerical integration to proceed to final flows.
The compressible model is intended to be the simplest possible model of a time-dependent magma delivery system, although no specific application is in mind. It omits variations in volatiles and viscosity, but it has these three important elements: 1. There is a single reservoir driven by a steady influx of material. The reservoir accumulates pressure to drive the melt upward through the colder surface of the earth. The reservoir in this model is linearly compressible, but that compressibility is meant to replace all the effects of buoyancy force driven by the density difference between magma and rock as well as the excess pressure from the elastic surroundings as magma accumulates under the region. 2. There is a permanent pathway to the surface, represented here by a simple cold tube with the added feature that it allows seepage flow. The pathway in our model represents a variety of natural pathways that guide magma ascent. There are cracks from stress in the elastic plate (abundantly observed seismically), the presence of brittle and weak material that develops cracks easily, and preheated aseismic pathways. 3. The melt can solidify along the tube. There are no volatiles, flow is one-dimensional low Re flow with composition and viscosity constant. Many magmatic systems (especially lavas) have laminar flow as used here (Dragoni et al., 2002;Klingelhofer et al., 1999;Rubin, 1993;Sakimoto & Gregg, 2001;Sakimoto & Zuber, 1998). Most important, the model eruption cannot happen unless the outflow is rapid enough to melt back the solid sheath of the tube (like the classic melt-back of a fissure as in Bruce & Huppert, 1989).
The dynamics of the spacing of active multiple tubes and the relation between spacing and the scaled manifold radius C −1/4 over a wide range of the values of C is obviously caused by the relatively close correlation between active flowing tube radius and scaled manifold radius. There is a small influence on spacing by the actual value of C. One might expect that these relations will not be the same for temperature-dependent viscosity laminar flow.
Although flow and freezing shut with true solidification differs from flow with viscosity variation, we found that invoking a minimum radius makes solidifying flows very similar to flow of fluids with large temperature-dependent viscosity. For example, when our minimum radius is inserted, there is a branch of the pressure curve that bends down to zero as flow approaches zero ( Figure 5 inset), just like flows with temperature-dependent viscosity (Helfrich, 1995;Whitehead & Helfrich, 1991;Wylie, Helfrich, et al., 1999;Wylie & Lister, 1995). Possibly the model of Wylie and Lister (1995) with a step change in viscosity is the closest equivalent to our solidification model, although their equations do not include a latent heat of fusion. In any case, the flow with variable viscosity inherently has a cold plug flow limit that is similar to seepage flow WHITEHEAD 10.1029/2020JF005969 19 of 25 so that our new results seem to apply to such problems. A systematic investigation of hysteresis might be more usefully conducted with viscosity variation.
A minimum radius in both Sections 3 and 5 is fundamental if one wants to avoid the discontinuity of freezing shut. For Section 3, the flow rate-pressure curve must have two extrema as in the inset of Figure 5. In that way, two intersections are stable and the third middle one is not, so oscillations can occur. In Section 5, the minimum radius prevents excessively large pressures that are associated with very small flow rates and vanishing radii.
Perhaps other physical processes such as temperature depending viscosity or temperature dependent terms in the heat equation can be invoked numerically instead of a minimum radius for solidifying flows in a future study. In any event, the need to invoke a minimum radius makes the results with large viscosity variation and with solidification very similar so future projects might simply use one or the other, depending on which is most convenient. In addition, some numerical results in Section 3 clearly apply to flow with viscosity variation and this should also be true for Section 5.
There are innumerable interesting extensions. One can combine these upstream conditions to flows with both viscosity variation and solidification, or have a slightly porous solid, or incorporate non-Newtonian flows like those reviewed by Kavanagh et al. (2018), or make a model of sedimentation problems or extend this approach to more complex flow geometry. It is not difficult to imagine the occurrence of very complicated or even truly chaotic flows. With enough complications, even realistic random-appearing patterns (Klein, 1982) could probably be generated. It is hoped, however, that the interesting behavior of these models with relatively simple flow situations can start to explain some of the elaborate piles of material that are encountered in igneous, frozen and depositional structures in the earth.

For small Stefan number (St
, the analysis of Holmes (2007) holds and there is some hope that those (and these) calculations might tie into natural situations. The values seem promising for many problems involving water, ice and magma. For water in a glacier,  / 0.0125 p H C L and a typical temperature difference between an intrusive liquid and the solidus of a few degrees has small St. This also might be true for some glacial drainage situations, although generally the value of Reynolds number and frictional heat generation would be great enough to suggest revisions of Equations A2 and A3. The revisions for turbulent transport of momentum and heat might alter the trend of pressure drop toward infinity at the approach of zero volume flux. The existing solutions (Holmes, 2007;Zerkle & Sunderland, 1968) seem to approach small Re in that limit. Although we can guess that flow should become laminar as glacial flows approach very low flux rate, a full study is still well-warranted. For magma,  / 0.0025 p H C L , and a magma temperature one or two hundred degrees above solidus has small St. This temperature difference is typical of most eruptions (values from Turcotte & Schubert, 2002). Volcanic eruptions producing hot water moving through tubular channels in ice do not fit the small St criterion.
What might these results imply for the spacing of outflows in nature? Let us try first to look at the formation of vents along a volcanic fissure. The number of tubes for N = 1,000 in Figure 10 is roughly fit by the relation Make a model of a fissure composed of 10 3 tubes spaced L m = 10 m apart feeding melt up from a shallow reservoir at a depth L = 1,000 m below the surface. A total flux of Q = 1 m 3 s −1 is evenly distributed at 1,000 m depth and therefore the flux per tube is Q = 0.001 m 3 s −1 . Using    2 / u n q Q LT along with magma thermal diffusivity  = 5 × 10 −7 m 2 s −1 , and T n = 10, gives u q = 0.020. As a first guess, equating the radius of the manifold tube to the tube going up to the surface so that r o = r m and C = L/L m then Equation 51 gives that # = 16/10 1/2 = 5 tubes that are active over the 10 km extent so there is a vent every 2 km. With greater depth of the fissure and everything else the same, then q u is smaller and there are fewer active tubes with wider vent spacing. These distances are plausible and given the great differences between this simple model and complex reality, the test seems to be promising.
Let us try a second example-the general problem of magma focusing at mid-ocean spreading centers. Pretend that there is a manifold consisting of a continuous mushy zone along a 1,000 km long ridge with vertical tubes each spaced 1 km apart that might bring melt up to the surface. To pick flux rate, we need to produce a flux that generates an oceanic crust thickness of 7 km with a ridge with a moderate spreading WHITEHEAD 10.1029/2020JF005969 20 of 25 rate of 0.1 m y −1 (=3.2 × 10 −9 m s −1 ). This gives a flux rate per tube spaced over the 1 km width covered by each tube of approximately 0.022 m 3 s −1 . Using a value of L = 30 km, (a minimum value for the depth), the same values of thermal diffusivity and T n as above, then the dimensionless value of flux rate is q u = 0.0149. There is little knowledge of what the equivalent for r m would be for either mushy zones or magma chambers under the ocean floor, so for a crude start use C = 1 (Note that a new model with a porous manifold is quite feasible.) This gives 12 active tubes for the ridge, equivalent to spacing of 83 km. This exceeds the spacing that is more typically 20-40 km for moderate rate mid-ocean spreading centers. Note also that this calculation implies that spacing is inversely proportional to flux rate so that with the present parameters ultraslow spreading centers might have spacings over 100 km and the fastest might have spacing less than 50 km.
These results help explain why magma cannot rise up everywhere in fissures and along spreading centers.
There are presumably ranges of parameters where volcanic intrusions might even freeze shut. Note that the volume flux rate used here is equal to 0.7 km 3 /y for the 1,000 km ridge, which reduces to a volume flux rate for each of the 12 tubes of 0.058 km 3 /y. This value is in the middle of the range of active volcanos in White et al. (2006) although they suggest other dynamics for governing the size of the volcanos. There are many other suggested dynamical factors governing the spacing of volcanos. To name a very few, there is Rayleigh-Taylor buoyancy that involves viscosity of the mushy zone (Schouten et al., 1985), there are combined buoyant, tectonic and mantle-forced flows (Magde & Sparks, 1997), and there is even deeper mantle flow (VanderBeek et al., 2016 and references therein). Results of this simple model suggest that lateral migration in the mushy zone with rising modulated by localized freezing dynamics might also be important and these dynamics can be added to the existing list.

Conclusions
A simple model is developed for liquid flowing into a freezing tube. This is used in conjunction with three different configurations of upstream flow to find parameters for instability leading to freezing.
The three configurations and their stability parameters are: • An upstream reservoir of finite size is fed steadily so pressure can change with upstream surface elevation. This has the dimensionless storage rate • Approaching very small values of each of these parameters produces freezing for an upstream flow less than three. Approaching very large values of each parameter has no freezing. • Numerical runs frequently have a difficulty when freezing is approached that is removed by allowing a minimum radius in the tube. Then, runs successfully occur over all time with all values of the parameters. • With multiple tubes, the active tube spacing is proportional to 1/4 C /q u with the others seeping. Therefore, an active tube adopts a resistance that is proportional to the cross manifold resistance. calculated profiles of   , I q . To calculate I, the solution for liquid temperature  uses the eigenvalues and eigenfunctions from Graetz (1883) as discussed further in Holmes (2007) and Holmes-Cerfon and Whitehead (2011). The results uniformly show that the radius of the solid surface decreases in the downstream distance (Figure 2c), and for steady flow, the scaling dictates that for given values of Q and  the dimensional distance downstream from the origin scales like  / q. This geometric independence from T n means the four interface profiles in Figure 2c are all versions of the same curve. The pressure drop across the tube is a function of volume flux rate q with a minimum value over the entire range, and p approaches infinity at the limits   0, q .