Numerical simulations of block-and-ash flows using the Titan2D flow model: examples from the 2006 eruption of Merapi Volcano,Java, Indonesia

We present Titan2D simulations of two well-characterized block-and-ash flow (BAF) events of the 2006 eruption of Merapi (Java, Indonesia) that affected the Gendol valley on the volcano’s southern flank and adjacent, densely populated interfluve (non-valley) areas: (1) a single dome-collapse event to the south that generated one of the smaller, post-June 14 flows and (2) a sustained, multiple dome-collapse event, also directed to the south, that produced the largest flows of the 2006 eruption emplaced in the afternoon of June 14. Using spatially varying bed friction angles, Titan2D is capable of reproducing the paths, velocities, runout distance, areas covered and deposited volumes of these flows over highly complex topography. The model results provide the basis for estimating the areas and levels of hazards associated with BAFs generated during relatively short as well as prolonged dome-collapse periods and guidance during future eruptive crises at Merapi.     

Field observations and surface characteristics of pristine block-and-ash flow deposits from the 2006 eruption of Merapi Volcano,Java, Indonesia

A series of pristine block-and-ash flow deposits from the May–June 2006 eruption of Merapi represent an exceptional record of small-volume pyroclastic flows generated by gravitational lava-dome collapses over a period of about two months. The deposits form nine overlapping lobes reaching ~ 7 km from the summit in the Gendol River valley on the volcano's southern flank, which were produced by successive flows generated during and after the major dome-collapse event on June 14. Both, single pulse (post-June 14 events) and multiple-pulse pyroclastic flows generated by sustained dome collapses on June 14 are recognised and three types of deposits, spread over an area of 4.7 km², are distinguished, totalling 13.3 × 10⁶ m³: (1) valley-confined basal avalanche deposits (11.7 × 10⁶ m³) in the Gendol River valley, (2) overbank pyroclastic-flow and associated surge deposits (1.4 × 10⁶ m³), where parts of the basal avalanche spread laterally onto interfluves and were subsequently channeled into the surrounding river valleys and (3) dilute ash-cloud surge deposits (0.2 × 10⁶ m³) along valley margins. Variations in the distribution, surface morphology and lithology of the deposits are related to the source materials involved in individual pyroclastic-flow-forming events and varying modes of transport and deposition of the different flows. Inferred flow velocities of the largest block-and-ash flows generated on June 14 vary from 43.8–13.5 m/s for the basal avalanche and from 62.6–24.2 m/s for the ash-cloud surge. The minimum temperatures range from 400 °C for the basal avalanche to 165 °C for the overlying ash cloud. Due to the potential of being re-channeled into adjacent river valleys and flowing laterally away from the main river channel, the overbank pyroclastic flows are considered the most hazardous part of the block-and-ash flow system. The conditions that lead to their development during flow transport and deposition must be taken into account when assessing future pyroclastic flow hazards at Merapi and similar volcanoes elsewhere.     

A three-dimensional Eulerian model for transport and deposition of volcanic ashes

Volcanic ash fallout represents a serious threat to people living near active volcanoes because it can produce several undesirable effects such as collapse of roofs by ash loading, respiratory sickness, air traffic disruption, or damage to agriculture. The assessment of such volcanic risk is therefore an issue of vital importance for public safety and its mitigation often requires to evaluate the temporal evolution of the phenomenon through reliable computational models.We develop an Eulerian model, named FALL3D, for the transport and deposition of volcanic ashes. The model is based on the advection–diffusion–sedimentation equation with a turbulent diffusion given by the gradient transport theory, a wind field obtained from a meteorological limited area model (LAM) and the source term derived from by buoyant plume theory. It can be used to forecast either ash concentration in the atmosphere or ash loading on the ground. Model inputs are topography, meteorological data given by a LAM, mass eruption rate, and a particle settling velocity distribution. A test application to the July 2001 Etna eruption is presented.     

A model based on Hirano-Exner equations for two-dimensional transient flows over heterogeneous erodible beds

In order to study the morphological evolution of river beds composed of heterogeneous material, the interaction among the different grain sizes must be taken into account. In this paper, these equations are combined with the two-dimensional shallow water equations to describe the flow field. The resulting system of equations can be solved in two ways: (i) in a coupled way, solving flow and sediment equations simultaneously at a given time-step or (ii) in an uncoupled manner by first solving the flow field and using the magnitudes obtained at each time-step to update the channel morphology (bed and surface composition). The coupled strategy is preferable when dealing with strong and quick interactions between the flow field, the bed evolution and the different particle sizes present on the bed surface. A number of numerical difficulties arise from solving the fully coupled system of equations. These problems are reduced by means of a weakly-coupled strategy to numerically estimate the wave celerities containing the information of the bed and the grain sizes present on the bed. Hence, a two-dimensional numerical scheme able to simulate in a self-stable way the unsteady morphological evolution of channels formed by cohesionless grain size mixtures is presented. The coupling technique is simplified without decreasing the number of waves involved in the numerical scheme but by simplifying their definitions. The numerical results are satisfactorily tested with synthetic cases and against experimental data.     

A 2D well-balanced,coupled model of water flow,sediment transport,and bed evolution based on unstructured grids with efficient variable storage strategy

In recent decades,a few Godunov-type,finite volume two-dimensional(2D)unstructured grid,coupled flow,and sediment models(GF2DUCM)have been developed for flows over erodible beds.These kinds of models are generally analyzed as a Vertex Model(VM)that define topography at the cell vertex,which can lead to the non-conservation of mass regarding flow,sediment,and bed evolution.Here,a full cellcantered variable storage method(Central Model or CM)is applied as the solution of the GF2DUCM.In this method,terrain elevation is defined at the cell centroids;this accurately describes the physical relations between the water depth and topography deformation.This approach can fully eliminate calculation errors in topography deformation at local cells caused by the interpolation of topography deformation at the cell vertex,and reduced uncertainty in the computation of the GF2DUCM.The model performance is systematically tested using a series of laboratory experiments,which demonstrate the mass conservation feature and high accuracy in reproducing hydrodynamic and morphological processes.     

A bed‐load transport model for rough turbulent open‐channel flows on plane beds

Data from flume studies are used to develop a model for predicting bed‐load transport rates in rough turbulent two‐dimensional open‐channel flows moving well sorted non‐cohesive sediments over plane mobile beds. The object is not to predict transport rates in natural channel flows but rather to provide a standard against which measured bed‐load transport rates influenced by factors such as bed forms, bed armouring, or limited sediment availability may be compared in order to assess the impact of these factors on bed‐load transport rates. The model is based on a revised version of Bagnold's basic energy equation i_bs_b = e_bω, where i_b is the immersed bed‐load transport rate, ω is flow power per unit area, e_b is the efficiency coefficient, and s_b is the stress coefficient defined as the ratio of the tangential bed shear stress caused by grain collisions and fluid drag to the immersed weight of the bed load. Expressions are developed for s_b and e_b in terms of G, a normalized measure of sediment transport stage, and these expressions are substituted into the revised energy equation to obtain the bed‐load transport equation i_b = ω G ^3·4. This equation applies regardless of the mode of bed‐load transport (i.e. saltation or sheet flow) and reduces to i_b = ω where G approaches 1 in the sheet‐flow regime. That i_b = ω does not mean that all the available power is dissipated in transporting the bed load. Rather, it reflects the fact that i_b is a transport rate that must be multiplied by s_b to become a work rate before it can be compared with ω. It follows that the proportion of ω that is dissipated in the transport of bed load is i_bs_b/ω, which is approximately 0·6 when i_b = ω. It is suggested that this remarkably high transport efficiency is achieved in sheet flow (1) because the ratio of grain‐to‐grain to grain‐to‐bed collisions increases with bed shear stress, and (2) because on average much more momentum is lost in a grain‐to‐bed collision than in a grain‐to‐grain one. Copyright © 2006 John Wiley & Sons, Ltd.     

A multiphase flow and multispecies reactive transport model for DNAPL-involved Compound Specific Isotope Analysis

This study presents a multiphase flow and multispecies reactive transport model for the simultaneous simulation of NAPL and groundwater flow, dissolution, and reactive transport with isotope fractionation, which can be used for better interpretation of NAPL-involved Compound Specific Isotope Analysis in 3D heterogeneous hydrogeologic systems. The model was verified for NAPL-aqueous phase equilibrium partitioning, aqueous phase multi-chain and multi-component reactive transport, and aqueous phase multi-component transport with isotope fractionation. Several illustrative examples are presented to investigate the effect of DNAPL spill rates, degradation rate constants, and enrichment factors on the temporal and spatial distribution of the isotope signatures of chlorinated aliphatic hydrocarbon groundwater plumes. The results clearly indicate that isotope signatures can be significantly different when considering multiphase flow within the source zone. A series of simulations indicate that degradation and isotope enrichment compete with dissolution to determine the isotope signatures in the source zone: isotope ratios remain the same as those of the source if dissolution dominates the reaction, while heavy isotopes are enriched in reactants along groundwater plume flow paths when degradation becomes dominant. It is also shown that NAPL composition can change from that of the injected source due to the partitioning of components between the aqueous and NAPL phases even when degradation is not allowed in NAPL phase. The three-dimensional simulation is presented to mechanistically illustrate the complexities in determining and interpreting the isotopic signatures with evolving DNAPL source architecture.     

An assembly model for simulation of large-scale ground water flow and transport

When managing large-scale ground water contamination problems, it is often necessary to model flow and transport using finely discretized domains--for instance (1) to simulate flow and transport near a contamination source area or in the area where a remediation technology is being implemented; (2) to account for small-scale heterogeneities; (3) to represent ground water-surface water interactions; or (4) some combination of these scenarios. A model with a large domain and fine-grid resolution will need extensive computing resources. In this work, a domain decomposition-based assembly model implemented in a parallel computing environment is developed, which will allow efficient simulation of large-scale ground water flow and transport problems using domain-wide grid refinement. The method employs common ground water flow (MODFLOW) and transport (RT3D) simulators, enabling the solution of almost all commonly encountered ground water flow and transport problems. The basic approach partitions a large model domain into any number of subdomains. Parallel processors are used to solve the model equations within each subdomain. Schwarz iteration is applied to match the flow solution at the subdomain boundaries. For the transport model, an extended numerical array is implemented to permit the exchange of dispersive and advective flux information across subdomain boundaries. The model is verified using a conventional single-domain model. Model simulations demonstrate that the proposed model operated in a parallel computing environment can result in considerable savings in computer run times (between 50% and 80%) compared with conventional modeling approaches and may be used to simulate grid discretizations that were formerly intractable.     

A model for distributed GIUH‐based flow routing on natural and anthropogenic hillslopes

Attempts to reduce the number of parameters in distributed rainfall–runoff models have not yet resulted in a model that is accurate for both natural and anthropogenic hillslopes. We take on the challenge by proposing a distributed model for overland flow and channel flow based on a combination of a linear response time distribution and the hillslope geomorphologic instantaneous unit hydrograph (GIUH), which can be calculated with only a digital elevation model and a map with field boundaries and channel network as input. The spatial domain is subdivided into representative elementary hillslopes (REHs) for each of which we define geometric and flow velocity parameters and compute the GIUH. The catchment GIUH is given by the sum of all REH responses. While most distributed models only perform well on natural hillslopes, the advantage of our approach is that it can also be applied to modified hillslopes with for example a rectangular drainage network and terrace cultivation. Tests show that the REH‐GIUH approach performs better than classical routing functions (exponential and gamma). Simulations of four virtual hillslopes suggest that peak flow at the catchment outlet is directly related to drainage density. By combining the distributed flow routing model with a lumped‐parameter infiltration model, we were also able to demonstrate that terrace cultivation delays the response time and reduces peak flow in comparison to the same hillslope, but with a natural stream network. The REH‐GIUH approach is a first step in the process of coupling distributed hydrological models to erosion and water quality models at the REH (associated with agricultural management) and at the catchment scale (associated with the evaluation of the environmental impact of human activities). It furthermore provides a basis for the development of models for large catchments and urban or peri‐urban catchments. Copyright © 2013 John Wiley & Sons, Ltd.     

Parametrisation of the orographic enhancement of precipitation and deposition in a long-term, long-range transport model

Orographic enhancement of wet deposition arising from the seeder-feeder effect is, by necessity, highly parametrised in long-range transport models of acid deposition that are long-term (i.e. annual average) and spatially resolved at tens of kilometres. Here, we describe a mechanistic approach to the incorporation of these mechanisms into such a model. The model formulation required the following: precipitation rate by direction and quantification of the fractions that are orographic and non-orographic; treatment of the fast oxidation of sulfur dioxide in clouds; the directionality of the seeder-feeder process; and a quantitative basis for increasing wet deposition factors to account for the seeder-feeder process. The directionality of non-orographic precipitation was determined from meteorological data at 47 sites across the UK. Orographic precipitation varies on a much finer scale than can be interpolated from measurements, and thus a modelling approach was adopted. The directionality of the seeder-feeder effect was taken from measurements. The enhancement factor of the orographic component of precipitation, assumed to represent feeder-rain, was determined from a review of measurements. Fast oxidation of sulfur dioxide is an observed phenomenon in cap-cloud, but limited in duration. An adjustment was made to the sulfur dioxide oxidation rate in the model in locations where cap-cloud was assumed to be present. The results from the model were compared with UK deposition budgets and enhanced wet deposition maps. The Revised parametrisation underestimated the UK wet deposition budgets of oxidised N and S, but spatial patterns of deposition were improved for much of the UK. It was concluded that this was a satisfactory outcome given the constraints of the statistical approach of weighting of deposition at receptors utilising straight line trajectories. The sensitivity of the model to directional constraints of seeder-feeder enhancement was tested and it was concluded that a fairly narrow constraint resulted in similar estimations to a broader one, and the broader constraint was thus adopted as frontal conditions which result in the process arrive from a fairly broad band of directions. When enhancement was allowed to occur from all directions, UK wet deposition of oxidised N and S was increased by 10%. The sensitivity to the enhancement factor on wet deposition was tested and found to be relatively robust. An increase in the enhancement factor from 2 to 6 resulted in increases in UK wet deposition of oxidised N and S of 9 and 6%, respectively.     

A parallel mesh-free contaminant transport model based on the Analytic Element and Streamline Methods

This paper introduces a new method for simulating large-scale subsurface contaminant transport that combines an Analytic Element Method (AEM) groundwater flow solution with a split-operator Streamline Method for modeling reactive transport. The key feature of the method is the manner in which the vertically integrated AEM flow solution is used to construct three-dimensional particle tracks that define the geometry of the Streamline Method. The inherently parallel nature of the algorithm supports the development of reactive transport models for spatial domains much larger than current grid-based methods. The applicability of the new approach is verified for cases with negligible transverse dispersion through comparisons to analytic solutions and existing numerical solutions, and parallel performance is demonstrated through a realistic test problem based on the regional-scale transport of agricultural contaminants from spatially distributed sources.     

A numerical model for the dynamics of pyroclastic flows at Galeras Volcano,Colombia

This paper presents a two-dimensional model for dilute pyroclastic flow dynamics that uses the compressible Navier–Stokes equation coupled with the Diffusion–Convection equation to take into account sedimentation. The model is applied to one of the slopes of Galeras Volcano to show: (1) the temperature evolution with the time; (2) dynamic pressure change; and (3) particle concentration along the computer domain from the eruption to the impact with a topographic barrier located more than 16 km from the source. Two initial solid volumetric fractions are modeled. For both cases, some of the structures located more distant than 10 km could survive, but in all cases the flow remains deadly. This paper shows that a dynamical model of pyroclastic flows can be implemented using personal computers.     

A semi-implicit three-dimensional numerical model for non-hydrostatic pressure free-surface flows on an unstructured,sigma grid

A semi-implicit 3-D numerical formulation for solving non-hydrostatic pressure free-surface flows on an unstructured,sigma grid is proposed.Pressure-splitting and 9 semi-implicit methods are inherited and reformed from Casulli’s z-coordinate formulation.The non-orthogonal sigma-coordinate transformation leads to additional terms. The resulting linear system for the non-hydrostatic correction is diagonally dominant but unsymmetric,and it is solved by the BiCGstab method.In contrast with z-coordinate non-hydrostatic models,the new model fits vertical boundaries much better,which is important for the long-time simulation of sediment transport and riverbed deformation.A lock-exchange density flow is computed to determine whether the new scheme is able to simulate non-hydrostatic free-surface flows.The new model is further verified using the field data of a natural river bend of the lower Yangtze River.Good agreement between simulations and earlier research results,field data is obtained, indicating that the new model is applicable to hydraulic projects in real rivers.     

A simple model accounting for the uptake,transport, and deposition of wind‐eroded mineral particles in the hyperarid coastal Atacama Desert of northern Chile

As previously observed in marine sediments collected downwind of African or South American continental sources, recent studies of sediment cores collected at the bottom of Mejillones Bay in north Chile (23°S) show a laminated structure in which the amount of particles of aeolian origin and their size create significant differences between the layers. This suggests inter‐annual to inter‐decadal variations in the strength of the local southerly winds responsible for (1) the erosion of the adjacent hyperarid surface of the Mejillones Pampa, and (2) the subsequent transport of the eroded particles towards the bay. A simple model accounting for the vertical uptake, transport, and deposition of the particles initially set into motion by wind at the surface of the pampa is proposed. This model, which could be adapted to other locations, assumes that the initial rate of (vertical) uptake is proportional to the (horizontal) saltation flux quantified by means of White's equation, that particles are lifted to a height (H), increasing with the magnitude of turbulence, and that sedimentation progressively removes the coarsest particles from the air column as it moves towards the bay. In this model, the proportionality constant (A) linking the vertical flux of particles with the horizontal flux, and the injection height (H) control the magnitude and size distribution of the deposition flux in the bay. Their values are determined using the wind speed measured over the pampa and the size distribution of particles collected in sediment traps deployed in the bay as constraints. After calibration, the model is used to assess the sensitivity of the deposition flux to the wind intensity variations. The possibility of performing such quantitative studies is necessary for interpreting precisely the variability of the aeolian material in the sediment cores collected at the bottom of Mejillones Bay. Copyright © 2010 John Wiley & Sons, Ltd.     

The simulation of sediment transport and erosion caused by free-surface flow based on two-phase SPH model with the improved Shields criterion

Ocean Dynamics - In this paper, a two-phase model based on the SPH method is established using the Newtonian fluid and non-Newtonian fluid to describe water and sediment respectively. In this...     

A new conceptual model for whole mantle convection and the origin of hotspot plumes

A new conceptual model of mantle convection is constructed for consideration of the origin of hotspot plumes, using recent evidence from seismology, high-pressure experiments, geodynamic modeling, geoid inversion studies, and post-glacial rebound analyses. This conceptual model delivers several key points. Firstly, some of the small-scale mantle upwellings observed as hotspots on the Earth's surface originate at the base of the mantle transition zone (MTZ), in which the Archean granitic continental material crust (TTG; tonalite-trondhjemite-granodiorite) with abundant radiogenic elements is accumulated. Secondly, the TTG crust and the subducted oceanic crust that have accumulated at the base of MTZ could act as thermal or mechanical insulators, leading to the formation of a hot and less viscous layer just beneath the MTZ; which may enhance the instability of plume generation at the base of the MTZ. Thirdly, the origin of some hotspot plumes is isolated from the large low shear-wave velocity provinces (LLSVPs) under Africa and the South Pacific. I consider that the conceptual model explains why almost all the hotspots around Africa are located above the margins of the African LLSVP. Because a planetary-scale trench system surrounding a “Pangean cell” has been spatially stable throughout the Phanerozoic, a large amount of the oceanic crustal layer is likely to be trapped in the MTZ under the Pangean cell. Therefore, under Africa, almost all of the hotspot plumes originate from the base of the MTZ, where a large amount of TTG and/or oceanic crusts has accumulated. This conceptual model may explain the fact that almost all the hotspots around Africa are located on margins above the African LLSVP. It is also considered that some of the hotspot plumes under the South Pacific thread through the TTG/oceanic crusts accumulated around the bottom of the MTZ, and some have their roots in the South Pacific LLSVP while others originate from the MTZ. The numerical simulations of mantle convection also speculate that the Earth's mantle convection is not thermally double-layered at the ringwoodite to perovskite + magnesiowüstite (Rw → Pv + Mw) phase boundary, because of its gentle negative Clapeyron slope. This is in contrast with some traditional images of mantle convection that have independent convection cells between the upper and lower mantle. These numerical studies speculate that the generation of stagnant slab at the base of the MTZ (as seismically observed globally) may not be due to the negative Clapeyron slope, and may instead be related to a viscosity increase (i.e., a viscosity jump) at the Rw → Pv + Mw phase boundary, or to a chemically stratified boundary between the upper and the lower mantle, as suggested by a recent high-pressure experiment.     

A conceptual hydrologic model for a forested Carolina bay depressional wetland on the Coastal Plain of South Carolina,USA

This paper describes how climate influences the hydrology of an ephemeral depressional wetland. Surface water and groundwater elevation data were collected for 7 years in a Coastal Plain watershed in South Carolina USA containing depressional wetlands, known as Carolina bays. Rainfall and temperature data were compared with water‐table well and piezometer data in and around one wetland. Using these data a conceptual model was created that describes the hydrology of the system under wet, dry, and drought conditions. The data suggest this wetland operates as a focal point for groundwater recharge under most climate conditions. During years of below‐normal to normal rainfall the hydraulic gradient indicated the potential for groundwater recharge from the depression, whereas during years of above‐normal rainfall, the hydraulic gradient between the adjacent upland, the wetland margin, and the wetland centre showed the potential for groundwater discharge into the wetland. Using high‐resolution water‐level measurements, this groundwater discharge condition was found to hold true even during individual rainfall events, especially under wet antecedent soil conditions. The dynamic nature of the hydrology in this Carolina bay clearly indicates it is not an isolated system as previously believed, and our groundwater data expand upon previous hydrologic investigations at similar sites which do not account for the role of groundwater in estimating the water budget of such systems. Copyright © 2007 John Wiley & Sons, Ltd.