Adaptive Modelling of Coupled Hydrological Processes with Applications in Water Economy

This project is structured into four sub projects. This group contributes with the design and implementation of a flexible and efficient software infrastructure for the solution of coupled systems of differential equations. Within the context of a domain decomposition approach this includes the coupling of equations, discretizations and linear solvers. The development is based on the software framework DUNE.

The prediction of flood waters requires the consideration of large areas and very complex local geometries. Hence, the parallel computability of the resulting problems is of great importance for this project. The storage and visualization of the computed results inevitably requires corresponding capabilities for data I/O and parallel rendering.

Dynamic Capillary Fringes

Modelling of Transport in strongly heterogeneous soils

Parameter Estimation for heterogeneous porous media

Development of an efficient and flexible Software Framework for partial differential equations (DUNE)

Most finite element, or finite volume software is built around a fixed mesh data structure. Therefore, each software package can only be used efficiently for a relatively narrow class of applications. For example, implementations supporting unstructured meshes allow the approximation of complex geometries but are in general much slower and require more memory than implementations using structured meshes.

By developing algorithms based on abstract interfaces, it is possible to avoid these restrictions, thus allowing multiple implementations of the same interface. Modern methods of generic programming (e.g. templates in C++) allow an efficient realization of this concept without restricting the flexibility of the resulting code. A new framework for partial differential equations - the "Distributed and Unified Numerics Environment" (DUNE) - is being developed according to this paradigm and in cooperation with many other universities. Currently, DUNE is applied in a vast variety applied sciences including the research on porous media, neuro transmitters, and particle accelerators.

The Iterative Solver Template Library

Sparse matrices obtained from finite element discretisations exhibit a lot of structure (e. g. discretisation of three-component system with linear finite elements and point-wise ordering) already known at compile time. As the knowledge is known already at compile time it can be exploited for efficiency using generic programming like in C++. The Iterative Solver Template Library (ISTL) provides a generic matrix/vector interface supporting a recursive block structure.

Solving large sparse linear systems is an ubiquitous task in the numerical solution of partial differential equations (PDEs). Increasing demands of computationally challenging applications both in problem size and algorithm complexity have lead to the development of parallel scalable solver libraries for these tasks. One of the most effective ways to achieve scalability is the use of multigrid or multilevel techniques. Algebraic Multigrid (AMG) is a very efficient algorithm for for solving large sparse problems on unstructured grids.

Our parallel AMG algorithm based on aggregation is designed to exploit the block structure of the matrices. Thus it can cope efficiently with scalar matrices as well as coupled block representations for systems of PDEs during the setup phase as well as during the solve phase. The algorithm proves to be a robust, efficient and scalable preconditioner within Krylow methods for the simulation of flow through heterogeneous media.

As a next step we want to tune the algorithm for solving linear systems stemming from Discontinuous Galerkin discretisations. Preliminary tests show that our algorithm is a potentially robust and efficient solver for these systems.

Unfitted Discontinuous Galerkin Methods

Simulation of physical, biological and chemical processes often involve complex shaped domains. Common problems are flow through root networks, solute transport on the pore scale of porous media or exchange processes through cell membranes.

Classical numerical methods require a grid resolving the complex geometry. Creating such grids is a very sophisticated process and therefore methods without this requirement are of great interest.

In this work a new approach to simulations on complex shaped domains was developed. The method is based on a Discontinuous Galerkin (DG) method with trial and test functions defined on a structured grid. Thus the number of degrees of freedom is proportional to the number of elements in the structured grid. The support of the trial and test functions is restricted according to the shape of the geometry. Essential Boundary conditions are imposed weakly via the Discontinuous Galerkin formulation. This method offers a discretization where the number of unknowns is independent of the complexity of the domain.

The method was successfully applied to stationary as well as time dependent problems.

Two Phase Navier Stokes Flow in Complex Domains

In this project, a robust and efficient method for the simulation of the laminar flow of two incompressible and immiscible fluids within a complex geometry is developed. The method for the solution of the corresponding Navier Stokes equations is based on the Unfitted Discontinuous Galerkin method which has been specifically developed for applications in complex domains. This nonconforming method does not require a global grid generation and realizes local conservation of impulse and mass. Furthermore, the pressure jump across the two phase interface can be approximated by the discontinuous base functions without introducing any numerical error due to the regularization of the pressure field near the interface.

The target applications of this model are the determination of hydraulic parameters for porous media multi phase flow models on a continuous macro scale. Especially the most recent models which take into account the immobile parts of the fluids and their interfacial area could greatly benefit from reliable simulations on the micro scale. The analysis of phenomena like hysteresis and fingering, which still defy all attempts of reliable modelling on the macro scale, are also part of this project.

Geometries of interesting soil samples are provided by the cooperation partner by means of x-ray tomography.

Large-scale numerical simulation of processes during CO2 storage in geologic formations

Geostatistical Inversion of Coupled Processes in Heterogeneous Porous Media

The assessment of the three-dimensional distribution of hydraulic parameters in groundwater bodies is of uttermost importance for the management of groundwater resources used for water supply, for the evaluation of anthropogenic impacts on water bodies including the loading with contaminants, and for the design of remediation schemes at contaminated sites. The distribution of hydraulic conductivity determines capture zones, flow paths, and travel times. Spatial variability occurs on practically all scales so that traditional conceptual models of the subsurface, in which the domain is subdivided into a few layers and zones with uniform coefficients, may be put into question. In geostatistical characterization, by contrast, the hydraulic parameters are considered as correlated random space variables. This allows applying inversion schemes, in which the estimated parameter fields continuously vary in space, meet the measurements of dependent quantities, and show the required spatial correlation. State-of-the-art geostatistical inversion schemes on serial computers can be used to estimate up to one million parameters. This, however, is hardly sufficient for a full three-dimensional representation of the heterogeneous subsurface. Additional restrictions of serial computing concern the inversion of large sets of measurements and the generation of multiple realizations in conditional Monte-Carlo simulations. The latter restrictions hamper the full use of geophysical data in joint inversion of hydraulic and geophysical surveys, and the use of steady-state concentration measurements for the estimation of hydraulic-conductivity fields.

The objective of the proposed project is to develop a program environment for geostatistical inversion of data originating from coupled flow, transport and geophysical assessment processes in heterogeneous porous formations using high-performance-computing techniques. As platform for the discretization and solution of forward and adjoint partial differential equations for groundwater flow, solute transport, and geophysical monitoring, the software framework "Distributed and Unified Numerics Environment" (DUNE), developed by the parallel computing group at University of Heidelberg and others will be adopted and improved. The geostatistical inversion methods developed by the hydrogeology group at University of Tübingen will be integrated into this framework, allowing fully parallelized computations for inversion and conditional Monte-Carlo simulations on three levels: domain decomposition, parallel evaluation of sensitivities, and parallel generation of conditional realizations.

The software framework developed in the proposed research will push the capabilities of jointly analyzing groundwater data from hydraulic experiments (heads, concentrations, flowmeter and pumping tests) and geophysical surveying (mainly geoelectrical tomography) onto a new level, and thus helps making better decisions for groundwater management and remediation under uncertainty.

Simulation of large-scale cortical network dynamics of morphologically detailed and electrophysiologically realistic neurons

Focus of this project is the simulation of networks on level of cortical units. A cortical unit of e.g. the rat barrel cortex consists of several thousand neurons. It is assumed that an individual column is involved in decision making processes as the gap crossing behavior of rodents.

In a first effort a simulator is realized that is able to generate and simulate neuron networks using a semi-statistical approach. The resulting diffusion-reaction type equation systems are solved in a scalable way using HPC resources. The forward simulator builds the basis to investigate the network dynamics inherent in individual realizations of semi-statistical generated networks. Furthermore inverse modelling will be important to fit measured electrophysiological network activity to simulated model network activity. Since no definitive information about the innervation is present we investigate several innervation scenarios in relation to the network response. Network activation is analyzed with regard to biophysically relevant quantities such as subthreshold response, number of activated synapses as well as spike output of the VPM-activated network.

The previously described statistical approach to build networks necessitates Monte-Carlo simulation to determine quantities of interest by ensemble averaging. An ensemble individual consists hereby of a single large-scale network realisation consisting of morphologically detailed neurons activated by the driving VPM cells with multivariate statistical wiring. To perform simulations and analyse responses of distinct ensemble series we use the numerical framework, NeuroDUNE, that has been developed to enable modeling and simulation of signal processing in such large-scale, full-compartmental neuron networks on sub cellular basis.

Reliable Simulation of Neurons and Networks

Signal processing in nerve cells can be modeled by non-linear systems of coupled ordinary and parabolic partial differential equations of reaction-diffusion type. Reliable simulations guarantee that errors in quantities of interest, such as point-wise approximation errors of the electrical potential, can be bounded by user-specified tolerances. Thus, reliable numerical solution schemes have to be able to detect spatial and temporal regions of high error contributions to arbitrarily chosen error functionals and to adapt spatial grids and time step sizes accordingly. In this project, reliable, robust and adaptive schemes for the simulation of single neurons and neuron networks are developed.

Controlling numerical errors in a robust manner require the use of mathematically sophisticated error estimators for both spatial and temporal error contributions. In a first step, error estimators developed in related research fields are adapted and employed in the scope of Computational Neuroscience. Their applicability and reliability is analyzed using various well known single cell nerve models. In a second step, the reliable single cell adaptive algorithms are employed in network simulations. As these networks are simulated using parallel computers, additional problems evolve from adaptivity. Amongst them are load balancing heuristics as well as guaranteed synchronicity of temporal events such as emission of action potentials.

All adaptive methods and error estimators described above have are realized using NeuroDUNE, a numerical simulation framework for the simulation of single neurons and large scale neuron networks. The adaptive, reliable algorithms are planned to become an integral part of NeuroDUNE.

Efficient and scalable 3D reconstruction methods for large-scale SBFSEM-datasets

In the last years novel techniques, such as the Serial Block-Face Scanning Electron Microscopy (SBFSEM) in medical and biological image acquisition allow to produce image data with a resolution in the range of several nanometers. Furthermore, the ongoing automatization of the image retrival process enables to scan considerable amounts of biological tissue in reasonable amounts of time and effort. Thus detailed, high quality image information is available and can be used to reconstruct nano-structures with physiological or anatomical relevance.

Focus of this work is the development of efficient and scalable algorithms for segmentation of SBFSEM image data, that allow a 3D reconstruction of individual neurons and finally of sample neural circuits. Emphasis is put on a robust and automatic reconstruction process that is able to work on large-scale data sets consisting of terabytes of information. Furthermore, this project focuses on the large-scale quantitative analysis, classification in types and topological cartography of spines in different dendritic branches, e.g. tuft vs. basal.

Main application area for the presented algorithms is data of mammial cortical structures. As example the extracted microanatomical information is relevant to decode simple decision making processes as the gap crossing behaviour of rodents.

Parameter estimation for realistic neuron channel models and tuning of large-scale network behavior

Efficient automatic parameter estimation techniques will be developed to (i) adapt single neuron channel kinetics, distributions, and voltage responses to match somatic and dendritic whole-cell recordings, and (ii) to adjust a limited set of network parameters such as synaptic connectivity strengths and radii to reproduce network level phenomena. For parameter estimation on the single-cell level we will employ constrained gradient-based methods with parameters ordered according to their sensitivities derived by numerical differentiation, a technique sufficient for a moderate number of parameters (< 30) which we have used successfully in other contexts. At the network level, gradient-free techniques such as simulated annealing and genetic algorithms will be used to estimate a limited set of moments of synaptic parameter distributions. To obtain an objective function at the network level, a number of crucial distributional and time series measures will be derived from the experimental data, such as interspike-interval and firing rate distributions.

Computational models for LFP signals

To connect biophysical models to EEG experiments, it is crucial to develop modules for LFP signal generation. This involves the coupling of network models to a model that describes electrical signal transduction through complicated three-dimensional extracellular space (ES). To represent the properties and complicated structure of extracellular space we build a model under usage of the DUNE framework. Two approaches are to be followed: (I) Modelling properties of ES transduction by equations with effective parameters. These parameters effectively take into account measured relations between multiple-unit spiking activity and simultaneously measured LFP and EEG signals. Using the DUNE framework will allow a direct coupling of the LFP model with the network simulations performed with NeuroDUNE. (II) A second method can consider the reconstructed ES per se. Such ES reconstructions have been performed by several researchers. Thus, a model motivated by first principles in complicated ES is build up from a pure mechanistic understanding of LFP. To take explicitly into account the distorted structure of ES specialized discretization and geometric modeling techniques are available within DUNE. Using experimental measurements the model will be refined and fine-tuned to allow predictions of LFP.