GIS and MATLAB Integrated for Groundwater Modelling

Bernard Raterman, Frans W. Schaars and Martin Griffioen

Abstract

This paper presents the applications of an integrated GIS and MATLAB system in several groundwater studies for Dutch Watercompanies. A recent study in the province of Utrecht proved that time and money can be saved by applying these integrated systems instead of using a conventional approach with existing groundwater modelling software. The integration of GIS and MATLAB is based on open dataformats and standards and may be used for many other application areas, especially where modelling in 3D and 4D is involved.


Introduction

Kiwa Research and Consultancy is the knowledge centre for water and associated ecological and environmental questions in The Netherlands. Kiwa is responsible for scheduling and conducting the joint research for Dutch water companies. Over the past 30 years expertise has been built up in diverse areas including hydrology, ecology, process technology and distribution technology. Since 1989 also GIS-technology has been used to support our core competence: the integration of these various fields of knowledge to develop innovative tools and concepts for adequate support of our customers.

In our daily work we have to deal with a broad range of data sources and formats and very different questions to solve. In the Netherlands a wide variety of groundwater modelling software is used. Among the finite-difference models the most commonly used is Modflow. Finite-element models Triwaco and MicroFem are also widespread. Manipulations are often necessary to exchange data from one model to the other. Therefore flexibility of our software tools is very important.

This paper illustrates why we have chosen to face some groundwater modelling challenges using GIS and MATLAB.

 

MATLAB

Some hydrologists have chosen not to be 'confined' to existing software for numerical groundwater modelling. They want to experiment more with analytical solutions for specific problems. Kees Maas and Theo Olsthoorn were pioneers in using MATLAB software (the MathWorks) to solve simple and complex hydrological problems in their daily work (Maas et al., 1997).

MATLAB stands for 'MATrix LABoratory'. It handles a range of computing tasks in engineering and science, from data acquisition and analysis to application development. The MATLAB environment integrates mathematical computing, visualization, and a powerful technical language. It is used in a variety of application areas, including signal and image processing, control system design, earth and life sciences, finance and economics, and instrumentation.

Because the data is stored in matrices of multiple dimensions a quick data access in 3D and 4D is also possible.

In many ways the evolution of MATLAB resembles the evolution of GIS, although in most application areas it has made a quicker transition from the universities to the industries. However, this does not hold true for applications in hydrology. This is due to the fact that hydrologists have built their models with the existing (numerical) modelling software. All their site-specific knowledge has been fed into these models, making them very valuable. There is an ongoing trend however to put this knowledge into GIS so they become available to a broader group of (potential) users.

Some problems may be tackled effectively with GIS techniques whereas others are more suitable for processing with MATLAB tools (e.g.3D and 4D modelling and visualization/animation). For these purposes a lot of AML's, Avenue scripts and MATLAB m-files have been developed over the past years.

 

The integration of ArcView and MATLAB

Integration levels

Heinzer et al. distinguishes three levels of integration ranging from model centric to GIS centric approaches. Partial integration uses GIS only for certain aspects of model parameter input. Encapsulated Integration is when GIS 'wraps' around a model and a GIS user and data-interface for the model exists. Embedded integration is when the solution algorithms are programmed within the GIS language.

We have used the embedded form in case of our hydro-ecological model NICHE® (Nature Impact of Changes in the Hydrological Environment) which predicts effects on vegetation after hydrological measures or changes in landuse. Actually this embedded integration was a logical step because no technical model and software implementation existed at the time we started developing it in 1994.

Throughout the years we have solved many little and a few big hydrological modelling problems with a project-specific partial or encapsulated integration using ArcInfo, ArcView and MATLAB. This was only possible after thorough discussions about which part of the processing or modelling could better be done either by GIS or by MATLAB. Only then it became clear which data to exchange, when and in which format. Although the results were satisfactory most of the times a general approach was still missing and we started to look for means to integrate on another level.

Towards a generic integration

We realized that we would be able to make a generic interface which we could effectively use in various (future) projects once we developed a simple method to define the 'links' we wanted, make m-files to read and write shapes and grids and establish simple means of communication between both applications. The close cooperation between GIS expert and hydrologist/MATLAB expert has led to the development of a valuable tool: the ArcView - MATLAB interface.

Components of the ArcView - MATLAB interface

The components currently making up the integration of ArcView and MATLAB are:

  1. A table describing all the necessary exchange parameters per MATLAB function (matlab_functions.dbf). Newly created functions are added to this table. In fact this is the simplest but most important core of the total integration.
  2. A self-adapting interface in ArcView produced with Dialog Designer and Avenue scripting. The dialog automatically changes with the selected function and corresponding definitions in "matlab_functions.dbf". Note: by using the "Switchyard programming technique" we were able to put all the Avenue code into one single script.
  3. The file av2mat.txt telling MATLAB what to do once the "RUN" option is chosen from the interface.
  4. The file mat2av.txt telling ArcView what to do after MATLAB has finished execution. This file is created by MATLAB.
  5. Two MATLAB m-files: "readshp.m" and "makeshp.m". These are conversion routines based on the open SHP standard. In MATLAB the data of the SHP, SHX and DBF files is kept in so-called 'structures'. GRIDS are currently exchanged as ASCII files and are stored as matrices in MATLAB.

A simple example

We will first use a simple example to illustrate how the integration works. Suppose we want to shift a point shape file with a certain dX and dY. First define what data we need to exchange; this is a point SHP and scalars for dX and dY. Therefore we add 3 records to the "matlab_functions.dbf" defining this new function which we call "Shift_point"

Added records in matlab_funtions.dbf

Figure 1: Records added to define the Shift_point function.

The [matfunc] field is used to call the appropriate MATLAB m-file, the [parameter] field is used to name the variables used in Avenue as well as in the m-files. The field [shape] tells ArcView and MATLAB what kind of data to expect (currently: scalar, grid, point, polyline, polygon and attributes of shapes). The field [descript] contains the description which will be displayed in the interface. Already the ArcView dialog will be able to display this function and when a shapefile of type point is available in the active view it will be listed in the combobox. The user can fill out the dX and dY values in the textline controls on the dialog.

Screendump of ArcView-MATLAB interface

Figure 2: The user interface with the active View displayed on top.

Now a MATLAB m-file called "Shift_point.m" is required doing the XY-shift when the "RUN" button is activated. Because it is not available yet we will create it:

Contents of Shift_point.m

file Shift_point.m

In this case the original shapefile will be overwritten. Naturally a new shapefile may also be produced instead. If the RUN option is chosen and all parameter values are defined the file av2mat.txt is generated to instruct MATLAB what to do.

Contents of av2mat.txt

file av2mat.txt

Now ArcView hands over the control to MATLAB. MATLAB interprets the above and carries out the job. When finished it writes mat2av.txt and the control is given back to ArcView.

Contents of mat2av.txt

file mat2av.txt

ArcView reads this line file. The word 'shp' is recognized as a keyword after which the complete pathname to the sourcefile is given. The result of the shift is shown in the next figure.

View with results

Figure 3: The results of the operation.

Other keywords may be 'grid' or 'error'. In the case of an error a message from MATLAB will display in an ArcView message box.

 

Case studies

Example 1) Water exploration in the province of Utrecht

Where can we withdraw good quality groundwater for the production of drinking water in the province of Utrecht without endangering groundwater-dependant vegetation? That was an important question in 1999 for Hydron Midden Nederland. This Dutch watercompany had to draw new plans for producing drinking water from new groundwater sources as some of the existing production wells were located in or near nature-conservation areas. Eventually these production wells will have to move due to government regulations. It was necessary to depict areas for further, more in-depth, studies within a few months time.

There was not enough time to run the existing numerical groundwater models on approximately 720.000 scenarios (eight capacity scenarios times 90.000 possible well locations). An analytical solution was already available in MATLAB. Since all the hydrological parameters of the confined aquifer system could easily be converted from the Triwaco models to polygon SHP-files, the most logical step was to combine these efforts into an integrated system ready to run this large number of scenarios. Together with the provincial government of Utrecht and Kiwa the following approach was defined.

  1. Indicate nature areas with groundwater dependant vegetation;
  2. Tolerate a maximum of 0.2 mm/day of seepage-loss in these areas;
  3. Define a maximum of 10 Million m3/year of groundwater to be withdrawn from the 3rd or 4th aquifer (approximately 60 to 160 m below surface level);
  4. 'Scan' the entire province for possible locations using the analytical solution in MATLAB based on a 200x200 meter grid (90.000 cells)

Figure 4 shows nature areas in the province of Utrecht and its vicinity. Three shades of green indicate the sensitivity of the vegetation to changes in groundwater table; dark green represents the most sensitive vegetation types.

Map of Utrecht and nature areas

Figure 4: The Province of Utrecht and selected nature areas.

The subsurface may be schematized as shown in figure 5. Here four waterbearing sandlayers are shown confined by clay and loam deposits. The yellow line demarcates the zone were the seepage loss becomes less then 0.2 mm per day. When the yellow line 'hits' the border of a nature area, the proposed well-site must be rejected. The distribution of all hydrological parameters were converted from polygon shapes (figure 6) into a 'stack' of GRID themes. These themes formed the input to the model in MATLAB. For each scenario a calculation was made using the ArcView - MATLAB integration.

Subsurface schematization

Figure 5: Schematization of the subsurface.

Thiessen polygons based on Triwaco model nodes

Figure 6: Hydrological parameters are stored as attributes of Thiessen polygons based on Triwaco model nodes.

After 16 hours of calculation time the final results were available in GIS. The map below shows a compilation of the eight output scenarios. In areas that are shaded dark-blue further studies can be undertaken to look for a site where 10 Million m3/year can be withdrawn without effecting nature areas. Lighter shades of blue indicate a potential annual capacity of 7.5, 5 or 2.5 Million m3.

Map with final results

Figure 7: Compilation of the results of all scenarios.

If conventional numerical modelling would have been used for this 'scan' the total processing time would have been much longer. Using the possibilities of ArcView and MATLAB we were able to produce a true 'quick-scan' and cut the costs of the project considerably.

Example 2) MATFLOW: MODFLOW in MATLAB via ArcView

The hydrologist accustomed to MATLAB will also be challenged to program MODFLOW in an m-file. Frans Schaars was and came up with 20 lines of code in Matflow.m. Input to MATFLOW must consist of regular grids. We have produced an example of a two layer model:

Definition of Matflow parameters

Figure 8: Parameter definition of the MATFLOW function.

We have used a polyline shapefile to model a river with a constant head (=3) and a point shape representing a well with a discharge of 1. Then we used Spatial Analyst to convert these shapes to regular grids and added some grids with constant values of Ones, Zeros, 0.1, 0.001 and also a Nodata-grid to serve as input to MATFLOW.

MATFLOW user interface

Figure 9: User interface to MATFLOW with the active View displayed on top.

The results after execution of MATFLOW with above parameters set is shown below.

Model output

Figure 10: The output of MATFLOW is put directly in ArcView.

We may use MATLAB again to look at a 3D surface plot of this drawdown pattern. This is a standard function in MATLAB and takes a grid as input.

MATLAB surfaceplot

Figure 11: Surfaceplot of the drawdown in layer 1 viewed from the northeast corner.

Example 3) SPREAD®

In 1998 the model SPREAD® (Solute Prediction from Regional distributed Agricultural Deposition) has been developed. The model predicts the mass-balance of nutrients (Nitrogen and Phosphorous) based on statistical data of manure and fertilizer use over the years combined with changes in land use, soil and hydrologic conditions. Concentrations of compounds in the soil and shallow-groundwater are predicted and may be used as input to transport models to estimate the quality of groundwater at an abstraction well. The model is implemented in ArcView and MATLAB. ArcView handles the data selection, polygon overlays, data aggregation and visualization. MATLAB routines predict the loads of the compounds per hectare, the concentrations in the unsaturated zone and in the phreatic groundwater taking into account biochemical processes such as denitrification.

Illustration of SPREAD

Figure 12: Illustration of the SPREAD®model.

Together with Artesia and RIZA (Institute for Inland Water Management and Waste Water Treatment, The Netherlands) we are currently working on extending SPREAD® to predict the surface water quality. Therefore interaction with other unsaturated zone, groundwater and surface water models is necessary.

 

Future plans

With the transition to ArcGIS 8.1 it becomes obvious that we want to rebuild the integration using the ActiveX technology. In doing so we expect to achieve even more flexibility. It will be possible to call ArcGIS COM objects from MATLAB and vice versa. The end-user will not have to be bothered with switching from application to application. We plan to extend the functionality of the interface e.g. dealing with multiple model-layers. We also want to extend the possibilities of the integration to facilitate 3D and 4D modelling applications.

 

Conclusions

For the time being the current methods used offer great possibilities, taking optimal advantage of both ArcView and MATLAB. Also the flexible interface tool built in Arcview is applicable to a wide variety of Avenue scripts developed over recent years. Once input (and output) are registered in a definition file like "matlab_functions.dbf" the interface may be used to address and run Avenue scripts.

The integration of GIS and MATLAB brings complex modelling closer to the end-user and the decision-maker. The hydrologist modelling in MATLAB benefits from presentation and accessiblility through ArcView while the GIS expert benefits from the powerful calculation, visualization and animation options offered by MATLAB.

 

References

Heinzer, T.J., M. Sebhat and B. Feinberg. 1999. Integrating Geographic Information Systems to Hydraulic Models: Methods and Case Studies. Esri User conference proceedings 1999.

Maas, C and Olsthoorn, T. 1997. Snelle oudjes gaan MATLAB (Golden oldies go MATLAB). Stromingen 3/4, p21-42.

N.V. Waterleidingbedrijf Midden-Nederland (Hydron Midden-Nederland). 1999. Vervangende Productie Capaciteit, bundeling technische deelrapporten. Technical reports 1 and 5.

The MathWorks, Inc. 2000. Using MATLAB, Version 6.


 

Author Information

Bernard Raterman
GIS expert, geologist
Kiwa Research and Consultancy
P.O. Box 1072
3430 BB Nieuwegein, The Netherlands
+31 30 6069541
Bernard.Raterman@kiwa.nl

 

Frans Schaars
Hydrologist
Artesia
Korte Weistraat 12
2871 BP Schoonhoven, The Netherlands
+31 182 387138
F.Schaars@artesia-water.nl

 

Martin Griffioen
Hydrologist
Hydron Midden-Nederland
P.O. Box 40205
3504 AA Utrecht, The Netherlands
+31 30 2487351
Mgriffioen@hydron-mn.nl