Past Mathematical Geoscience Seminar

The collapse of granular columns in a viscous fluid is a common model case for submarine geophysical flows. In immersed granular collapses, dense packings result in slow dynamics and short runout distances, while loose packings are associated with fast dynamics and long runout distances. However, the underlying mechanisms of the triggering and runout, particularly regarding the complex fluid-particle interactions at the pore-scale, are yet to be fully understood. In this study, a three-dimensional approach coupling the Lattice Boltzmann Method and the Discrete Element Method is adopted to investigate the influence of packing density on the collapsing dynamics. The direct numerical simulation of fluid-particle interactions provides evidence of the pore pressure feedback mechanism. In dense cases, a strong arborescent contact force network can form to prevent particles from sliding, resulting in a creeping failure behavior. In contrast, the granular phase is liquefied substantially in loose cases, leading to a rapid and catastrophic failure. Furthermore, hydroplaning can take place in loose cases due to the fast-moving surge front, which reduces the frictional resistance dramatically and thereby results in a longer runout distance. More quantitatively, we are able to linearly correlate the normalized runout distance and the densimetric Froude number across a wide range of length scales, including small-scale numerical/experimental data and large-scale field data.

Arctic sea ice is one of the most sensitive components of the Earth’s climate system. The underlying ocean plays an important role in the evolution of the ice cover through its heat flux at the ice-ocean interface. Despite its importance, the spatio-temporal variations of this heat flux are not well understood. In this talk, I will take the following approach to study the variations in the heat flux. First, I will consider the problem of classical Rayleigh-Bénard convection and systematically explore the effects of fractal boundaries on heat transport using direct numerical simulations. And second, I will analyze time-series data from the Surface Heat Budget of the Arctic Ocean (SHEBA) program using Multifractal Detrended Fluctuation Analysis (MFDFA) to understand the nature of fluctuations in the heat flux. I will also discuss developing simple stochastic ODEs using results from these studies.

Of the canonical stratified shear flow instabilities (Kelvin–Helmholtz, Holmboe-wave and Taylor–Caulfield), the Taylor–Caulfield instability (TCI) has received relatively little attention, and forms the focus of the presentation. A diagnostic of the linear instability dynamics is developed that exploits the net pseudomomentum to distinguish TCI from the other two instabilities for any given flow profile. Next, the nonlinear dynamics of TCI is shown across its range of unstable horizontal wavenumbers and bulk Richardson numbers. At small bulk Richardson numbers, a cascade of billow structures of sequentially smaller size may form. For large bulk Richardson numbers, the primary nonlinear travelling waves formed by the linear instability break down via a small-scale, Kelvin– Helmholtz-like roll-up mechanism with an associated large amount of mixing. In all cases, secondary parasitic nonlinear Holmboe waves appear at late times for high Prandtl number. Finally, a nonlinear diagnostic is proposed to distinguish between the saturated states of the three canonical instabilities based on their distinctive density–streamfunction and generalised vorticity–streamfunction relations.

We study the stability properties of the Eulerian mean flow generated by monochromatic surface-gravity waves propagating in a rotating frame. The wave averaged equations, also known as the Craik-Leibovich equations, govern the evolution of the mean flow. For propagating waves in a rotating frame these equations admit a steady depth-dependent base flow sometimes called the Ekman-Stokes spiral, because of its resemblance to the standard Ekman spiral. This base flow profile is controlled by two non-dimensional numbers, the Ekman number Ek and the Rossby number Ro. We show that this steady laminar velocity profile is linearly unstable above a critical Rossby number Roc(Ek). We determine the threshold Rossby number as a function of Ek using a numerical eigenvalue solver, before confirming the numerical results through asymptotic expansions in the large/low Ek limit. These were also confirmed by nonlinear simulations of the Craik-Leibovich equations. When the system is well above the linear instability threshold, Ro >> Roc, the resulting flow fluctuates chaotically. We will discuss the possible implications in an oceanographic context, as well as for laboratory experiments.

This work aims at developing new approach for parameterizing mesoscale eddy effects for use in non-eddy-resolving ocean circulation models. These effects are often modelled as some diffusion process or a stochastic forcing, and the proposed approach is implicitly related to the latter category. The idea is to approximate transient eddy flux divergence in a simple way, to find its actual dynamical footprints by solving a simplified but dynamically relevant problem, and to relate the ensemble of footprints to the large-scale flow properties.

The Kelvin wave is perhaps the most important of the equatorially trapped waves in the terrestrial atmosphere and ocean, and plays a role in various phenomena such as tropical convection and El Nino. Theoretically, it can be understood from the linear dynamics of a stratified fluid on an equatorial beta plane, which, with simple assumptions about the disturbance structure, leads to wavelike solutions propagating along the equator, with exponential decay in latitude. However, when the simplest possible background flow is added (with uniform latitudinal shear), the Kelvin wave (but not the other equatorial waves) becomes unstable. This happens in an extremely unusual way: there is instability for arbitrarily small nondimensional shear p, and the growth rate is proportional to exp(-1/p^2) as p->0. This in contrast to most hydrodynamic instabilities, in which the growth rate typically scales as a positive power of p-p_c as the control parameter p passes through a critical value p_c.

This Kelvin wave instability has been established numerically by Natarov and Boyd, who also speculated as to the underlying mathematical cause. Here we show how the growth rate and full spatial structure of the instability may be derived using matched asymptotic expansions applied to the (linear) equations of motion. This involves an adventure with Whittaker functions in the exponentially-decaying tails of the Kelvin waves, and a trick to reveal the exponentially small growth rate from a formulation that only uses regular perturbation expansions. Numerical verification of the analysis is also interesting and challenging, since special high-precision solutions of the governing ODE are required even when the nondimensional shear is not that small (circa 0.5).

In this talk, we will consider how two very different atmospheric phenomena, the terrestrial tropical cyclone and the martian polar vortex, can be described within a single simplified dynamical framework based on the forced shallow water equations. Dynamical forcings include angular momentum transport by secondary (transverse) circulations and local heating due to latent heat release. The forcings act in very different ways in the two systems but in both cases lead to distinct annular distributions of potential vorticity, with a local vorticity maximum at a finite radius surrounding a central minimum.  In both systems, the resulting vorticity distributions are subject to shear instability and the degree of eddy growth versus annular persistence can be examined explicitly under different forcing scenarios.

Inertia-gravity waves (IGWs) are ubiquitous in the ocean and the atmosphere. Once generated (by tides, topography, convection and other processes), they propagate and scatter in the large-scale, geostrophically-balanced background flow. I will discuss models of this scattering which represent the background flow as a random field with known statistics. Without assumption of spatial scale separation between waves and flow, the scattering is described by a kinetic equation involving a scattering cross section determined by the energy spectrum of the flow. In the limit of small-scale waves, this equation reduces to a diffusion equation in wavenumber space. This predicts, in particular, IGW energy spectra scaling as k^{-2}, consistent with observations in the atmosphere and ocean, lending some support to recent claims that (sub)mesoscale spectra can be attributed to almost linear IGWs.  The theoretical predictions are checked against numerical simulations of the three-dimensional Boussinesq equations.
(Joint work with Miles Savva and Hossein Kafiabad.)

Energy transfers from small-scale turbulence and waves to large-scale flows are ubiquituous in oceans, atmospheres, planetary cores and stars.

Therefore, turbulence and waves have a direct effect on the large-scale organization of geophysical and astrophysical fluids and can affect their long-term dynamics.

In this talk I will discuss recent direct numerical simulation (DNS) results of two upscale energy transfer mechanisms that emerge from the dynamics of a fluid that is self-organized in a turbulent layer next to a stably-stratified one. This self-organization in an adjacent "two-layer" turbulent-stratified system is ubiquituous in nature and is representative of e.g. Earth's troposphere-stratosphere system, the oceans' surface mixed layer-thermocline system, and stars' convective-radiative interiors. The first set of DNS results will demonstrate how turbulent motions can generate internal waves, which then force a slowly-reversing large-scale flow, akin to Earth's Quasi-Biennial Oscillation (QBO). The second set of DNS results will show how the stratified layer regulates the emergence of large-scale vortices (LSV) in the turbulent layer under rapid rotation in the regime known as geostrophic turbulence. I will demonstrate why it is important to resolve both the turbulence and the waves, as otherwise the natural variability of the QBO is lost and LSV cannot form. I will discuss future works and highlight how the results may guide the implementation of upscale energy transfers in global earth system models.

Partially molten materials resist shearing and compaction. This resistance

is described by a fourth-rank effective viscosity tensor. When the tensor

is isotropic, two scalars determine the resistance: an effective shear and

an effective bulk viscosity. In this seminar, calculations are presented of

the effective viscosity tensor during diffusion creep for a 3D tessellation of

tetrakaidecahedrons (truncated octahedrons). The geometry of the melt is

determined by assuming textural equilibrium.  Two parameters

control the effect of melt on the viscosity tensor: the porosity and the

dihedral angle. Calculations for both Nabarro-Herring (volume diffusion)

and Coble (surface diffusion) creep are presented. For Nabarro-Herring

creep the bulk viscosity becomes singular as the porosity vanishes. This

singularity is logarithmic, a weaker singularity than typically assumed in

geodynamic models. The presence of a small amount of melt (0.1% porosity)

causes the effective shear viscosity to approximately halve. For Coble creep,

previous modelling work has argued that a very small amount of melt may

lead to a substantial, factor of 5, drop in the shear viscosity. Here, a

much smaller, factor of 1.4, drop is obtained.

Ice streams are fast flowing regions of ice that generally slide over a layer of unconsolidated, water-saturated subglacial sediment known as till.  A striking feature that has been observed geophysically is that subglacial till has been found to accumulate distinctively into sedimentary wedges at the grounding zones (regions where ice sheets begin to detach from the bedrock to form freely floating ice shelves) of both past and present-day ice sheets. These grounding-zone wedges have important implications for ice-sheet stability against grounding zone retreat in response to rising sea levels, and their origins have remained a long-standing question. Using a combination of mathematical modelling, a series of laboratory experiments, field data and numerical simulations, we develop a fluid-mechanical model that explains the mechanism of the formation of these sedimentary wedges in terms of the loading and unloading of deformable till in three dynamical regimes. We also undertake a series of analogue laboratory experiments, which reveal that a similar wedge of underlying fluid accumulates spontaneously in experimental grounding zones, we formulate local conditions relating wedge slopes in each of the scenarios and compare them to available geophysical radargram data at the well lubricated, fast-flowing Whillans Ice Stream.

Many scientists, and in particular mathematicians, report difficulty in understanding thermodynamics. So why is thermodynamics so difficult? To attempt an answer, we begin by looking at the components in an exposition of a scientific theory. These include a mathematical core, a motivation for the choice of variables and equations, some historical remarks, some examples and a discussion of how variables, parameters, and functions (such as equations of state) can be inferred from experiments. There are other components too, such as an account of how a theory relates to other theories in the subject.

It will be suggested that theories of thermodynamics are hard to understand because (i) many expositions appear to argue from the particular to the general (ii) there are several different thermodynamic theories that have no obvious logical or mathematical equivalence (iii) each theory really is subtle and requires intense study (iv) in most expositions different theories are mixed up, and the different components of a scientific exposition are also mixed up. So, by presenting one theory at a time, and by making clear which component is being discussed, we might reduce the difficulty in understanding any individual thermodynamic theory. The key is perhaps separation of the mathematical core from the physical motivation. It is also useful to realise that a motivation is not generally the same as a proof, and that no theory is actually true.

By way of illustration we will attempt expositions of two of the simplest thermodynamic theories – reversible and then irreversible thermodynamics of homogeneous materials – where the mathematical core and the motivation are discussed separately. In conclusion we’ll relate these two simple theories to other, foundational and generalised, thermodynamic theories.

The accurate modelling of geophysical flows often requires information which is difficult to measure and therefore poorly quantified. Such information may relate to the fluid properties or an unknown boundary condition, for example. The premise of this talk is that when the flow is bounded by a free surface, the deformation of this free surface contains useful information which can be used to infer such unknown quantities. The increasing availability of free surface data through remote sensing using drones and satellites provides the impetus to develop new mathematical methods and numerical tools to interpret the signature embedded in the free surface deformation. This talk will explore two recent examples drawn from glaciology and inspired from volcanology for which free surface data was successfully used to reconstruct an unknown field.

The magma chamber - an underground vat of fluid magma that is tapped during volcanic eruptions - has been the foundation of models of volcanic eruptions for many decades and successfully explains many geological observations.  However, geophysics has failed to image the postulated large magma chambers, and the chemistry and ages of crystals in erupted magmas indicate a more complicated history.  New conceptual models depict subsurface magmatic systems as dominantly uneruptible crystalline networks with interstitial melt (mushes) extending deep into the Earth's crust to the mantle, containing lenses of potentially eruptible (low-crystallinity) magma.  These lenses would commonly be less dense than the overlying mush and so Rayleigh Taylor instabilities should develop leading to ascent of blobs of magma unless the growth rate is sufficiently slow that other processes (e.g. solidification) dominate.  The viscosity contrast between a buoyant layer and mush is typically extremely large; a consequence is that the horizontal dimension of a magma reservoir is commonly much less than the theoretical fastest growing wavelength assuming an infinite horizontal layer.  

I will present laboratory experiments and linear stability analysis for low Reynolds number, laterally confined Rayleigh Taylor instabilities involving one layer that is much thinner and much less viscous than the other.  I will then apply the results to magmatic systems, comparing timescales for development of the instability and the volumes of packets of rising melt generated, with the frequencies and sizes of volcanic eruptions.  I will then discuss limitations of this work and outstanding fluid dynamical problems in this new paradigm of trans-crustal magma mush systems.

Glacier flow is an example of a gravity driven non-linear viscous flow at low Reynolds numbers. As a glacier flows over an undulating bed, the surface topography is modified in response. Some information about bed conditions is therefore contained in the shape of the surface and the surface velocity field. I will present theoretical and numerical work on how basal conditions on glaciers affect ice flow, and how one can obtain information about basal conditions through surface-to-bed inversion. I’ll give an overview over inverse methodology currently used in glaciology, and how satellite data is now routinely used to invert for bed properties of the Greenland and the Antarctic Ice Sheets.

Plumes are a characteristic feature of convective flow through porous media. Their dynamics are an important part of numerous geological processes, ranging from mixing in magma chambers to the convective dissolution of sequestered carbon dioxide. In this talk, I will discuss models for the spread of convective plumes in a heterogeneous porous environment. I will focus particularly on the effect of thin, roughly horizontal, low-permeability barriers to flow, which provide a generic form of heterogeneity in geological settings, and are a particularly widespread feature of sedimentary formations. With the aid of high-resolution numerical simulations, I will explore how a plume spreads and flows in the presence of one or more of these layers, and will briefly consider the implications of these findings in physical settings.

In the Boussinesq framework, velocity couples to density fluctuations whereas in magnetohydrodynamic turbulence, the velocity field is coupled to the magnetic field. Both systems support waves (inertia-gravity in the presence of rotation, or Alfvén), with anisotropic dispersion relations. What kind of turbulence regimes result from the interactions between waves and nonlinear eddies in such flows? And what is delimiting these regimes?

I shall sketch the phenomenological framework for rotating stratified turbulence within which one is led to scaling laws in terms of the Froude number, Fr=U/[LN], which measures the relative celerity of gravity waves and nonlinear eddies, with U and L characteristic velocity and length scale, and N the Brunt-V\"ais\"al\"a frequency. These laws apply to the mixing efficiency of such flows, indicating the relative roles of the buoyancy flux due to the waves, and of the measured kinetic and potential energy dissipation rates. Various measures of mixing are found to follow power laws in terms of the Froude number, and may differ for the three regimes that can be identified, namely the wave-dominated, wave-eddy balance and eddy-dominated domains [1]. In particular, in the intermediate regime, the effective dissipation varies linearly with Fr, in agreement with simple wave-turbulence arguments. This analysis is inspired by and corroborates results from a large parametric study using direct numerical simulations (DNS) on grids of 1024^3 points, as well as from atmospheric and oceanic observations.

Such scaling laws can be related to previous DNS results concerning the existence for the energy of bi-directional constant-flux cascades to both the small scales and to the large scales, due to the presence of rotation in such flows, as measured for example in the ocean. These dual energy cascades lead to an alteration, and a decrease, of the mixing and available energy to be dissipated in the small scales [2]. Some perspectives might also be given at the end of the talk.

[1] A. Pouquet, D. Rosenberg, R. Marino & C. Herbert, Scaling laws for mixing and dissipation in unforced rotating stratified turbulence. J. Fluid Mechanics 844, 519, 2018.
[2] R. Marino, A. Pouquet & D. Rosenberg, Resolving the paradox of oceanic large-scale balance and small-scale mixing. Physical Review Letters 114, 114504, 2015.

A granular material forms a distinct and fascinating phase in physics -- sand acts as a fluid as grains flow through your fingers, the fallen grains form a solid heap on the floor or may suspend in the wind like a gas.

The main challenge of studying granular materials is the development of constitutive models valid across scales, from the micro-scale (collisions between individual particles), via the meso-scale (flow structures inside avalanches) to the macro-scale (dunes, heaps, chute flows).

In this talk, I am highlighting three recent projects from my laboratory, each highlighting physical behavior at a different scale. First, using the property of birefringence, we are quantifying both kinetic and dynamic properties in an avalanche of macroscopic particles and measure rheological properties. Secondly, we explore avalanches on an erodible bed that display an intriguing dynamic intermittency between regimes. Lastly, we take a closer look at aqueous (water-driven) dunes in a novel rotating experiment and resolve an outstanding scaling controversy between migration velocity and dune dimension.

A Surstseyan eruption is a particular kind of volcanic eruption which involves the bulk interaction of water and hot magma. Surtsey Island was born during such an eruption process in the 1940s. I will talk about mathematical modelling of the flashing of water to steam inside a hot erupted lava ball called a Surtseyan bomb. The overall motivation is to understand what determines whether such a bomb will fragment or just quietly fizzle out...
Partial differential equations model transient changes in temperature and pressure in Surtseyan ejecta. We have used a highly simplified approach to the temperature behaviour, to separate temperature from pressure. The resulting pressure diffusion equation was solved numerically and asymptotically to derive a single parametric condition for rupture of ejecta. We found that provided the permeability of the magma ball is relatively large, steam escapes rapidly enough to relieve the high pressure developed at the flashing front, so that rupture does not occur. This rupture criterion is consistent with existing field estimates of the permeability of intact Surtseyan bombs, fizzlers that have survived.
I describe an improvement of this model that allows for the fact that pressure and temperature are in fact coupled, and that the process is not adiabatic. A more systematic reduction of the resulting coupled nonlinear partial differential equations that arise from mass, momentum and energy conservation is described. We adapt an energy equation presented in G.K. Batchelor's book {\em An Introduction to Fluid Dynamics} that allows for pressure-work. This is work in progress.  Work done with Emma Greenbank, Ian Schipper and Andrew Fowler 

Through laboratory experiments, we examine the transport, settling and resuspension of sediments as well as the influence of floating particles upon damping wave motion.   Salt water is shown to enhance flocculation of clay and hence increase their settling rate.   In studies modelling sediment-bearing (hypopycnal) river plumes, experiments show that the particles that eventually settle through uniform-density fluid toward a sloping bottom form a turbidity current.  Meanwhile, even though the removal of particles should increase the buoyancy and hence speed of the surface current, in reality the surface current stops.  This reveals that the removal of fresh water carried by the viscous boundary layers surrounding the settling particles drains the current even when their concentration by volume is less than 5%. The microscopic effect of boundary layer transport by particles upon the large scale evolution is dramatically evident in the circumstance of a mesopycnal particle-bearing current that advances along the interface of a two-layer fluid.  As the fresh water rises and particles fall, the current itself stops and reverses direction.  As a final example, the periodic separation and consolidation of particles floating on a surface perturbed by surface waves is shown to damp faster than exponentially to attain a finite-time arrest as a result of efficiently damped flows through interstitial spaces between particles - a phenomenon that may be important for understanding the damping of surface waves by sea ice in the Arctic Ocean (and which is well-known to anyone drinking a pint with a proper head or a margarita with rocks or slush).

In this talk, I show how concepts from non-equilibrium statistical physics can be employed in the study of climate. The specific problem addressed is the geophysical-scale evolution of Arctic sea ice. Using an analogy with Brownian motion, the original evolution equation for the sea ice thickness distribution function by Thorndike et al. (J. Geophys. Res. 80(33), pp. 4501 — 4513, 1975) is transformed to a Fokker-Planck-like conservation law. The steady solution is $g(h) = {\cal N}(q) h^q \mathrm{e}^{-~ h/H}$, where $q$ and $H$ are expressible in terms of moments over the transition probabilities between thickness categories. The solution exhibits the functional form used in observational fits and shows that for $h \ll 1$, $g(h)$ is controlled by both thermodynamics and mechanics, whereas for $h \gg 1$ only mechanics controls $g(h)$. We also derive the underlying Langevin equation governing the dynamics of the ice thickness $h$, from which we predict the observed $g(h)$. Further, seasonality is introduced by using the Eisenman-Wettlaufer model (Proc. Natl. Acad. Sci. USA 106, pp. 28-32, 2009) for the thermal growth of sea ice. The time-dependent problem is studied by numerically integrating the Fokker-Planck equation. The results obtained from these numerical integrations and their comparison with satellite observations are discussed.

Numerical simulation provides an important contribution to the management of oil reservoirs, and the ‘reservoir simulator’ has been an essential tool for reservoir engineers since the 1970’s. I will describe the role of the ‘well model’ in reservoir simulation. Its main purpose is to determine the production and injection flows of the reservoir fluids at the surface under a variety of operating constraints, and to supply source and sink terms to the grid cells of the reservoir model.

Advances in well technology (horizontal, multilateral, and smart wells containing flow control devices) have imposed additional demands on the well model. It must allow the fluid mixture properties to vary with position in the well, and enable different fluid streams to comingle. Friction may make an important contribution to the local pressure gradient. To provide an improved representation of the physics of fluid flow, the well is discretised into a network of segments, where each segment has its own set of variables describing the multiphase flow conditions. Individual segments can be configured to represent flow control devices, accessing lookup tables or built-in correlations to determine the pressure drop across the device as a function of the flow conditions.

The ability to couple the wells to a production facility model such as a pipeline network is a crucial advantage for field development and optimization studies, particularly for offshore fields. I will conclude by comparing two techniques for combining a network model with the reservoir simulation. One method is to extend the simulator’s well model to include the network, providing a fully integrated reservoir/well/network simulation. The other method is to run the reservoir and facility models as separate simulations coupled by a ‘controller’, which periodically balances them by exchanging boundary conditions. The latter approach allows the engineer to use a choice of specialist facility simulators.