Search

The behaviour of low-viscosity, pressure-driven compressible pore fluid flows in viscously deformable porous media is studied here with specific application to gas flow in lava domes. The combined flow of gas and lava is shown to be governed by a two-equation set of nonlinear mixed hyperbolic–parabolic type partial differential equations describing the evolution of gas pore pressure and lava porosity. Steady state solution of this system is achieved when the gas pore pressure is magmastatic and the porosity profile accommodates the magmastatic pressure condition by increased compaction of the medium with depth. A one-dimensional (vertical) numerical linear stability analysis (LSA) is presented here. As a consequence of the pore-fluid compressibility and the presence of gravitation compaction, the gradients present in the steady-state solution cause variable coefficients in the linearized equations which generate instability in the LSA despite the diffusion-like and dissipative terms in the original system. The onset of this instability is shown to be strongly controlled by the thickness of the flow and the maximum porosity, itself a function of the mass flow rate of gas. Numerical solutions of the fully nonlinear system are also presented and exhibit nonlinear wave propagation features such as shock formation. As applied to gas flow within lava domes, the details of this dynamics help explain observations of cyclic lava dome extrusion and explosion episodes. Because the instability is stronger in thicker flows, the continued extrusion and thickening of a lava dome constitutes an increasing likelihood of instability onset, pressure wave growth and ultimately explosion.

Episodes of slow uplift and subsidence of the ground, called bradyseism, characterize the recent dynamics of the Campi Flegrei caldera (Italy). In the last decades two major bradyseismic crises occ...

Uncertainty in predictions from a model of volcanic ash transport in the atmosphere arises from uncertainty in both eruption source parameters and the model wind field. In a previous contribution, we analyzed the probability of ash cloud presence using weighted samples of volcanic ash transport and dispersal model runs and a reanalysis wind field to propagate uncertainty in eruption source parameters alone. In this contribution, the probabilistic modeling is extended by using ensemble forecast wind fields as well as uncertain source parameters. The impact on ash transport of variability in wind fields due to unresolved scales of motion as well as model physics uncertainty is also explored. We have therefore generated a weighted, probabilistic forecast of volcanic ash transport with only a priori information, exploring uncertainty in both the wind field and the volcanic source.

In this paper, we present new ideas to greatly enhance the quality of uncertainty quantification in the DDDAS framework. We build on ongoing work in large scale transport of geophysical mass of volcanic origin – a danger to both land based installations and airborne vehicles. The principal new idea introduced is the concept of a localized Bayes linear model as a surrogate for the expensive simulator. Probability of ash presence is compared to earlier work.

In this paper, we will present ongoing work on using a dynamic data driven application system (DDDAS) based approach to the forecast of volcanic ash transport and dispersal. Our primary modeling tool will be a new code puffin formed by the combination of a plume eruption model Bent and the ash transport model Puff. Data from satellite imagery, observation of vent parameters and windfields will drive our simulations. We will use ensemble based uncertainty quantification and parameter estimation methodology – polynomial chaos quadrature in combination with data integration to complete the DDDAS loop.

Current shallow granular flow models suited to arbitrary topography can be divided into two types, those formulated in bed-fitted curvilinear coordinates, and those formulated in global Cartesian coordinates. The shallow granular flow model of Denlinger and Iverson \cite{Denlinger2004} and the Boussinesq-type shallow granular flow theory of Castro-Orgaz \emph{et al}. \cite{Castro2014} are formulated in a Cartesian coordinate system (with $z$ vertical), and both account for the effect of nonzero vertical acceleration on depth-averaged momentum fluxes and stress states. In this paper, we first reformulate the vertical normal stress of Castro-Orgaz \emph{et al}. \cite{Castro2014} in a quadratic polynomial in the relative elevation $\eta$. This form allows for analytical depth integration of the vertical normal stress. We then calculate the basal normal stress based on the basal friction law and scaling analysis. These calculations, plus certain constitutive relations, lead to a refined full non-hydrostatic shallow granular flow model, which is further rewritten in a form of Boussinesq-type water wave equations for future numerical studies. In the present numerical study, we apply the open-source code TITAN2D to numerical solution of a low-order version of the full model involving only a mean vertical acceleration correction term. To cure the numerical instability associated with discretization of the enhanced gravity, we propose an approximate formula for the enhanced gravity by utilizing the hydrostatic pressure assumption in the bed normal direction. Numerical calculations are conducted for several test cases involving steep slopes. Comparison with a bed-fitted model shows that even the simplified non-hydrostatic Cartesian model can be used to simulate shallow granular flows over arbitrary topography., Comment: 29 pages

This paper describes a new methodology to quantify the variation in the output of a computational fluid dynamics model for block and ash flows, when the digital elevation model (DEM) of the terrain and other inputs are given as a range of possible values with a prescribed uncertainty. Integrating these variations in the possible flows as a function of input uncertainties provides well-defined hazard probabilities at specific locations, i.e. a hazard map. Earlier work provided a methodology for assessing hazards based on variations in flow initiation and friction parameters. This paper extends this approach to include the effect of terrain error and uncertainty. The results are based on potential flows at Mammoth Mountain, CA, and Galeras Volcano, Colombia. The analysis establishes the soundness of the approach and the effect of including the uncertainty in DEMs in the construction of probabilistic hazard maps.

In this paper we will present early work on using a DDDAS based approach to the construction of probabilistic estimates of volcanic ash transport and dispersal. Our primary modeling tools will be a combination of a plume eruption model BENT and the ash transport model PUFF. Data from satellite imagery, observation of vent parameters and windfields will drive our simulations. We will use uncertainty quantification methodology – polynomial chaos quadrature in combination with data integration to complete the DDDAS loop.

Many model-based investigation techniques, such as sensitivity analysis, optimization, and statistical inference, require a large number of model evaluations to be performed at different input and/or parameter values. This limits the application of these techniques to models that can be implemented in computationally efficient computer codes. Emulators, by providing efficient interpolation between outputs of deterministic simulation models, can considerably extend the field of applicability of such computationally demanding techniques. So far, the dominant techniques for developing emulators have been priors in the form of Gaussian stochastic processes (GASP) that were conditioned with a design data set of inputs and corresponding model outputs. In the context of dynamic models, this approach has two essential disadvantages: (i) these emulators do not consider our knowledge of the structure of the model, and (ii) they run into numerical difficulties if there are a large number of closely spaced input points as is often the case in the time dimension of dynamic models. To address both of these problems, a new concept of developing emulators for dynamic models is proposed. This concept is based on a prior that combines a simplified linear state space model of the temporal evolution of the dynamic model with Gaussian stochastic processes for the innovation terms as functions of model parameters and/or inputs. These innovation terms are intended to correct the error of the linear model at each output step. Conditioning this prior to the design data set is done by Kalman smoothing. This leads to an efficient emulator that, due to the consideration of our knowledge about dominant mechanisms built into the simulation model, can be expected to outperform purely statistical emulators at least in cases in which the design data set is small. The feasibility and potential difficulties of the proposed approach are demonstrated by the application to a simple hydrological model.

Debris flows, avalanches, landslides, and other geophysical mass flows can contain O(106–1010) m3 or more of material. These flows commonly consist of mixture of soil and rocks with a significant quantity of interstitial fluid. They can be tens of meters deep, and their runouts can extend many kilometers. The complicated rheology of such a mixture challenges every constitutive model that can reasonably be applied; the range of length and timescales involved in such mass flows challenges the computational capabilities of existing systems.This paper extends recent efforts to develop a depth averaged "thin layer" model for geophysical mass flows that contain a mixture of solid material and fluid. Concepts from the engineering community are integrated with phenomenological findings in geo-science, resulting in a theory that accounts for the principal solid and fluid forces as well as interactions between the phases, across a wide range of solid volume fraction. A principal contribution here is to present drag and phase interaction terms that comport with the literature in geo-sciences. The program predicts the evolution of the concentration and dynamic pressure. The theory is validated with with data from one dimensional dam break solutions and it is verified with data from artificial channel experiments.

This paper describes the development of highly accurate adaptive discontinuous Galerkin schemes for the solution of the equations arising from a thin layer type model of debris flows. Such flows have wide applicability in the analysis of avalanches induced by many natural calamities, e.g. volcanoes, earthquakes, etc. These schemes are coupled with special parallel solution methodologies to produce a simulation tool capable of very high-order numerical accuracy. The methodology successfully replicates cold rock avalanches at Mount Rainier, Washington and hot volcanic particulate flows at Colima Volcano, Mexico.

The collapse and decrepitation of a lava dome at the summit of a volcano generally results in the generation of dense granular flows, often referred to as block and ash flows. As the dome particles propagate from the source, they break apart by internal pressure as well as collision. The propagation of block and ash flows can be simulated to some accuracy with a depth averaged numerical model of the equations of continuity and momentum for a material with a frictional resistance. However, important features of such flows, such as the influence of remote stress through force chains, erosion of the volcano substrate, and shock formation and pressurization upon particle break up are poorly understood. In the near future, the influence of these factors will be incorporated into depth averaged models. Various numerical techniques based on particles will some day yield results that can be compared not only with bulk flow properties, but to the internal layering of block and ash flow deposits.

The Titan2D geophysical mass-flow model is evaluated by comparing its simulation results and those obtained from another flow model, FLOW3D, with published data on the 1963 Little Tahoma Peak avalanches on Mount Rainier, Washington. The avalanches, totaling approximately 10� 10 6 m 3 of broken lava blocks and other debris, traveled 6.8 km horizontally and fell 1.8 km vertically (H/L=0.246). Velocities calculated from runup range from 24 to 42 m/s and may have been as high as 130 m/s while the avalanches passed over Emmons Glacier. Titan2D is a code for an incompressible Coulomb continuum; it is a depth-averaged, dshallow-waterT, granular-flow model. The conservation equations for mass and momentum are solved with a Coulomb-type friction term at the basal interface. The governing equations are solved on multiple processors using a parallel, adaptive mesh, Godunov scheme. Adaptive gridding dynamically concentrates computing power in regions of special interest; mesh refinement and coarsening key on the perimeter of the moving avalanche. The model flow initiates as a pile defined as an ellipsoid by a height (z) and an elliptical base defined by radii in the x and y planes. Flow parameters are the internal friction angle and bed friction angle. Results from the model are similar in terms of velocity history, lateral spreading, location of runup areas, and final distribution of the Little Tahoma Peak deposit. The avalanches passed over the Emmons Glacier along their upper flow paths, but lower in the valley they traversed stream gravels and glacial outwash deposits. This presents difficulty in assigning an appropriate bed friction angle for the entire deposit. Incorporation of variable bed friction angles into the model using GIS will help to resolve this issue. D 2004 Elsevier B.V. All rights reserved.

High fidelity computational simulation can be an invaluable tool in planning strategies for hazard risk mitigation. The accuracy and reliability of the predictions are crucial elements of these tools being successful. We present here a new simulation tool for dry granular avalanches using several new techniques to create a highly accurate tool with a combination of good models and numerical methodologies. Highlights of our methodology are the use of a depth averaged model of the conservation laws and an adaptive grid Godunov solver to solve the resulting equations The software is designed to run on distributed memory supercomputers and make use of digital elevation data dynamically i.e. refine the grid and input data to finer resolutions to better capture flow features as the flow evolves. Our simulations are validated using quantitative and qualitative comparisons to tabletop experiments and data from field observations. Our software is freely available and can be used on a range of hardware from inexpensive desktop computers to the largest supercomputers.

We present a framework for modeling a dry geophysical mass of granular material -- a debris or volcanic avalanche or landslide -- flowing over an erodible surface. We also describe a computing environment that incorporates topographical data into a parallel, adaptive mesh computational algorithm that solves the model equations.

Method-dependent mechanisms that may affect dynamic numerical solutions of a hyperbolic partial differential equation that models concentration profiles in renal tubules are described. Some numerical methods that have been applied to the equation are summarized, and ways by which the methods may misrepresent true solutions are analysed. Comparison of these methods demonstrates the need for thoughtful application of computational mathematics when simulating complicated time-dependent phenomena.

Recent micropuncture studies in rats have demonstrated the existence of oscillatory states in nephron filtration mediated by tubuloglomerular feedback (TGF). We develop a minimal mathematical model of the TGF system, consisting of a first-order hyperbolic partial differential equation describing thick ascending limb (TAL) NaCl reabsorption and an empirical feedback relation. An analytic bifurcation analysis of this model provides fundamental insight into how oscillatory states depend on the physiological parameters of the model. In the special case of no solute backleak in the TAL, the emergence of oscillations explicitly depends on two nondimensional parameters. The first corresponds to the delay time of the TGF response across the juxtaglomerular apparatus, and the second corresponds to the product of the slope of the TGF response curve at the steady-state operating point and the space derivative of the steady-state NaCl concentration profile in the TAL at the macula densa. Numerical calculations for the case without TAL backleak are consistent with this result. Numerical simulation of the more general case with TAL backleak shows that the bifurcation analysis still provides useful predictions concerning nephron dynamics. With typical parameter values, the analysis predicts that the TGF system will be in oscillatory state. However, the system is near enough to the boundary of the nonoscillatory region so that small changes in parameter values could result in nonoscillatory behavior.

[1] This paper presents several standard and new methods for characterizing the effect of input data uncertainty on model output for hazardous geophysical mass flows. Note that we do not attempt here to characterize the inherent randomness of such flow events. We focus here on the problem of characterizing uncertainty in model output due to lack of knowledge of such input for a particular event. Methods applied include classical Monte Carlo and Latin hypercube sampling and more recent stochastic collocation, polynomial chaos, spectral projection and a newly developed extension thereof named polynomial chaos quadrature. The simple and robust samplings based Monte Carlo type methods are usually computationally intractable for reasonable physical models, while the more sophisticated and computationally efficient polynomial chaos method often breaks down for complex models. The spectral projection and polynomial chaos quadrature methods discussed here produce results of quality comparable to the polynomial chaos type methods while preserving the simplicity and robustness of the Monte Carlo-type sampling based approaches at much lower cost. The computational efficiency, however, degrades with increasing numbers of random variables. A procedure for converting the output uncertainty characterization into a map showing the probability of a hazard threshold being exceeded is also presented. The uncertainty quantification procedures are applied first in simple settings to illustrate the procedure and then subsequently applied to the 1991 block-and-ash flows at Colima Volcano, Mexico.

The tubuloglomerular feedback (TGF) system regulates the fluid and solute load in each nephron of the mammalian kidney. In this report, we obtain a necessary condition for sustained oscillations in a simple mathematical model for the TGF loop. This model consists of a hyperbolic partial differential equation representing the chloride concentration in the thick ascending limb and an empirical function describing the feedback response of the juxtaglomerular apparatus. A bifurcation analysis shows that critical parameters are the time delay, the slope of the empirical feedback relation at the steady state, and the steady-state space derivative of the chloride concentration at the macula densa.

The ability to dynamically change data input to a computation is a key feature enabling simulation to be used in many applications. In this study, computation of geophysical mass flow is updated on the fly by changing terrain data. Accommodating such changes in a parallel environment entails new developments in parallel data management and gridding. Adaptivity, and in particular unrefinement, is critical for maintaining parallel efficiency. The application under study in this work is the result of a multidisciplinary collaboration between engineers, mathematicians, geologists, and hazard assessment personnel. In addition, adaptive gridding enables efficient use of computational resources, allowing for run-time determination of optimal computing resources. Combining these attributes allows run time conditions to inform calculations, which in turn provide up-to-date information to hazard management personnel.

A transport-reaction model describing penetration of plasmin by diffusion and permeation into a dissolving fibrin gel was solved numerically to explore mechanisms that lead to the formation and growth of dissolution fingers through blood clots during thrombolytic therapy. Under conditions of fluid permeation driven by arterial pressures, small random spatial variations in the initial fibrin density within clots (+/-4 to 25% peak variations) were predicted by the simulation to result in dramatic dissolution fingers that grew in time. With in vitro experiments, video microscopy revealed that the shape of the proximal face of a fibrin gel, when deformed by pressure-driven permeation, led to lytic breakthrough in the center of the clot, consistent with model predictions of increased velocities in this region leading to cannulation. Computer simulation of lysis of fibrin retracted by platelets (where more permeable regions are expected in the middle of the clot due to retraction) predicted cannulation of the clot during thrombolysis. This residual, annular thrombus was predicted to lyse more slowly, because radial pressure gradients to drive inner clot permeation were quite small. In conjunction with kinetic models of systemic pharmacodynamics and plasminogen activation biochemistry, a two-dimensional transport-reaction model can facilitate the prediction of the time and causes of clot cannulation, poor reperfusion, and embolism during thrombolysis.

The load of water and solute entering each nephron of the mammalian kidney is regulated by the tubuloglomerular feedback (TGF) mechanism, a negative feedback loop. Experiments in rats have shown that key variables of this feedback system may exhibit TGF-mediated oscillations. Mathematical modeling studies have shown that the open-feedback-loop gain is a crucial parameter for determining whether oscillations will emerge. However, two different formulations of this gain have been used. The first is the steady-state gain, a readily measurable quantity corresponding to the steady-state reduction in single-nephron glomerular filtration rate (SNGFR) subsequent to a sustained increased in ascending limb flow rate. The second is an instantaneous gain, a variable arising from theoretical considerations corresponding to the maximum reduction in SNGFR resulting from an instantaneous shift of the ascending limb flow column, with the assumption that the SNGFR response is also instantaneous. Here we show by an analytic argument how the steady-state and instantaneous open-feedback-loop gains for the ascending limb are related. In the case of no solute backleak into the ascending limb, the two formulations of gain are equivalent; however, in the presence of solute backleak, the instantaneous gain is larger in magnitude than the steady-state gain. With typical physiological parameters for the rat, calculations with a model previously devised by us show that the gains differ by 5-10%. Hence, experimental measurements of the steady-state gain may provide useful lower-bound estimates of the instantaneous gain of the feedback system in the normal rat. However, the gains may diverge significantly in pathophysiological states where ascending limb transport is compromised by abnormally high NaCl permeability.

This project combines analytic and computational investigations to understand the dynamics of elastic-plastic deformation of granular materials, particularly issues related to the formation of shearbands. Roughly speaking, shearbands form when the governing equations cease to be wellposed. Our research examines the issue of well-posedness, loss of hyperbolicity, and regularization. This final report summarizes our work on (1) computation of deformation and formation of shear bands in granular material; (1) analysis of a gradient theory of granular plasticity; (2) related elastic and viscoelastic systems of PDE which may lose hyperbolicity.... Granular material, Plastic deformation, Hyperbolic equations.