Table of Contents

#Table-of-Contents
Loading output library...

THE STUDY AREA

#THE-STUDY-AREA

Leaflet slippy map

#Leaflet-slippy-map
Loading output library...

Set map extents

#Set-map-extents

Initial Cartopy plot of the study area

#Initial-Cartopy-plot-of-the-study-area
Loading output library...

Helper functions to construct a gradually-refined grid

#Helper-functions-to-construct-a-gradually-refined-grid
  • 3 levels of refinement (far-field, near-field, local-field)

STRATEGY:

  • Locate center point of innermost grid with the greater refinement (local-area)
  • Starting from the center point, add cells outwards until reaching the bounding box of the refinement area
  • gradually refine using 1.5 rule up to the following refinement step
  • repeat for upper levels of refinement
Loading output library...

Plot grid on leaflet map (uses mplleaflet)

#Plot-grid-on-leaflet-map-(uses-mplleaflet)

Create workspace directory, set path to MODFLOW executable, model name, etc.

#Create-workspace-directory,-set-path-to-MODFLOW-executable,-model-name,-etc.

Define model name of the model and the location of MODFLOW executable. All MODFLOW files and output will be stored in the subdirectory defined by the workspace. Create a model named ml and specify that this is a MODFLOW-2005 model.

STEADY STATE

#STEADY-STATE

Temporal discretisation

#Temporal-discretisation

Initial spatial discretisation

#Initial-spatial-discretisation

Create DIS package object

#Create-DIS-package-object

TRANSIENT

#TRANSIENT

Temporal discretisation

#Temporal-discretisation

Predevelopment period

t=0-1970: burn-in, steady state, development of saltwater wedge

t=1970-1980: dynamic cyclic equilibrium (no pumping)

Development period

t=1980-1990: pumping begins (no calibration)

t=1990-2016(Dec): calibration

Loading output library...
Loading output library...

LAYER ELEVATIONS

#LAYER-ELEVATIONS

Set spatial reference for model grid

#Set-spatial-reference-for-model-grid
Loading output library...

NOTE: The BAS object can be easily and intuitively constructed with GeoPandas. First we read in the MODFLOW grid shapefile created with Flop and create a GeoPandas DataFrame from it

Loading output library...

Pickle the geodataframe

Plot the grid to make sure everything is OK

Loading output library...

Model top elevations (DEM)

#Model-top-elevations-(DEM)
Loading output library...
Loading output library...
Loading output library...

Plot the DEM and histogram of elevations

Loading output library...
Loading output library...
Loading output library...
Loading output library...

DEMraster_to_MODFLOW helper function to map the DEM to MODFLOW grid

Check how the MODFLOW_gpd looks like at this stage

Loading output library...

Pickle geodataframe

Plot DEM raster as applied to MODFLOW grid

Loading output library...
Loading output library...

Now we need to overwrite the initial DIS package, using the new elevations set for model top

Double-check that L1_top was correctly created

Loading output library...
Loading output library...

Superficial aquifer — bottom elevations

#Superficial-aquifer-—-bottom-elevations
Loading output library...
Loading output library...
Loading output library...

Plot Superficial aquifer bottom elevations and historgram of elevations

Loading output library...

DEMraster_to_MODFLOW helper function to map the DEM to MODFLOW grid

Apply the function

Loading output library...

Plot DEM raster applied to MODFLOW grid

Loading output library...
Loading output library...

Now we need to overwrite the initial DIS package, using the new elevations for model top

Double-check that model_bottom was correctly created

Loading output library...
Loading output library...

Create final DIS package object

#Create-final-DIS-package-object
Loading output library...
Loading output library...
Loading output library...
Loading output library...

For clarity, reorganise the order of columns in the MODFLOW geodataframe

#For-clarity,-reorganise-the-order-of-columns-in-the-MODFLOW-geodataframe
Loading output library...

Initialise by setting all cells to 'inactive' --> ibound=0:

Loading output library...

Read and create GeoDataFrame for BCs

Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...

Spatial joins are used to map boundary conditions (shapefiles) to the MODFLOW grid

Add ibound values to MODFLOW geodataframe

Loading output library...
Loading output library...
Loading output library...

make sure rows and columns are ints for indexing purposes

Loading output library...

Pickle geodataframe

Loading output library...

Pickle geodataframe

create ibound array

Create BAS package object

Loading output library...
Loading output library...
Loading output library...
Loading output library...

Create array with stress period data for the GHB package

Loading output library...
Loading output library...
Loading output library...

Hydraulic conductivity zones

#Hydraulic-conductivity-zones
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...

The intersection operation between the shapefiles and MODFLOW grid may assign cells to two hydrogeological zones. We can arbitrarily resolve ambiguities like so:

Loading output library...

Pickle geodataframe

Homogeneous HK (one zone)

#Homogeneous-HK-(one-zone)

Create LPF package object

Homogeneous HK (multiple zones)

#Homogeneous-HK-(multiple-zones)
Loading output library...

Heterogeneous HK (pilot points + Radial Basis Functions)

#Heterogeneous-HK-(pilot-points-+-Radial-Basis-Functions)
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...

Heterogeneous HK (pilot points + KNN regressor)

#Heterogeneous-HK-(pilot-points-+-KNN-regressor)
Loading output library...
Loading output library...

STEADY STATE

#STEADY-STATE

Set parameters

#Set-parameters

Create EVT package object

#Create-EVT-package-object

TRANSIENT

#TRANSIENT

Read SILO historical records

#Read-SILO-historical-records

Display complete record

#Display-complete-record
Loading output library...

pre-1970 record

#pre-1970-record
Loading output library...

use pre-1970 mean evaporation for the burn-in period

Loading output library...
Loading output library...
Loading output library...

Make nper (stress period) the index

Loading output library...
Loading output library...

Change units (to L/T)

#Change-units-(to-L/T)

The maxET values are in mm/month -> need to convert to L/T (i.e., m/d; x 12/(365*1000)

Loading output library...

Create EVT package object

#Create-EVT-package-object

read drainage network shapefile; set correct CRS

Loading output library...

plot the shapefile

Loading output library...
Loading output library...

we need to do a spatial join with the MODFLOW grid geodataframe i.e., tag those cells that will perform as drains

Plot intersection drainage network mapped to MODFLOW grid

Loading output library...
Loading output library...

Recall main MODFLOW geodataframe

Loading output library...

Incorporate the drainage network to the geodataframe; checking that it was properly appended

Loading output library...

need to write stress period data for the DRN package; create a geodataframe with the table entires that intersect a drain

Loading output library...

Create the array with stress period data for the DRN package

Loading output library...

Create DRN package object

Loading output library...
Loading output library...

STEADY STATE

#STEADY-STATE

Create RCH package object

#Create-RCH-package-object

TRANSIENT

#TRANSIENT

Write MODFLOW input files and run model

#Write-MODFLOW-input-files-and-run-model
Loading output library...
Loading output library...
Loading output library...

Retrieve heads from .hds file

Fixed masked data (-999.98... to 999)

Plot heads and DTWT using ModelMap

#Plot-heads-and-DTWT-using-ModelMap
Loading output library...
Loading output library...

mask inactive cells

compute depth to water table (DTWT)

Loading output library...

Calibration will proceed using the Spotpy package:

http://fb09-pasig.umwelt.uni-giessen.de/spotpy/

A typical Spotpy calibration setup is as follows:

Calibration targets (GWL observations)

#Calibration-targets-(GWL-observations)

Load pickle

#Load-pickle
Loading output library...

Drop duplicates

#Drop-duplicates
Loading output library...

Quality check drilling depth of calibration target bores

#Quality-check-drilling-depth-of-calibration-target-bores
Loading output library...
Loading output library...

Drop records with drilling depth likely deeper than superficial aquifer

#Drop-records-with-drilling-depth-likely-deeper-than-superficial-aquifer
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...

WTCHM FORECASTING SCENARIOS

  • CLIMATE: min, max, ave (x2, CSIRO/DWER) = 6 scenarios
  • PUMPING RATES: current, additional (60%, 80%...of MAR flux) = 3 scenarios
  • PUMPING LOCATIONS: = 2 scenarios?
  • MAR RATES: = 3 scenarios

PERFORMANCE METRICS

  • SEAWATER INTRUSION
  • INUNDATION
  • NUTRIENT LOADS TO COCKBURN SOUND
  • TRAVEL TIMES (OCEAN, WETLANDS)
  • MAR WATER DILUTION (0-1)

Compute drain length in cell

#Compute-drain-length-in-cell

roadLengthInCell = roadsMultiLineString.intersection(cellRect).length

Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...
Loading output library...

Land Use

#Land-Use

Performs land use reclassification from MRS datasets

IMPROVEMENTS:

  • Parse year from filename
  • Save as pickles raster_year
  • return year (to pass on to plot function)
Loading output library...

Performs BATCH land use reclassification from MRS datasets

IMPROVEMENTS:

  • Parse year from filename
  • Save as pickles raster_year
  • return year (to pass on to plot function)

Tweak previous function to take year as an argument

Set pivot table path

Plot function to compare original and reclassified rasters

IMPROVEMENTS:

  • Check that raster and MODFLOW grid have the same CRS (GeoPandas complains if not the same, anyway...)
  • take year as argument
  • save as figures raster_year
Loading output library...

Applies raster data to a MODFLOW grid (via GeoPandas)

IMPROVEMENTS:

  • Check that raster and MODFLOW grid have the same CRS (GeoPandas complains if not the same, anyway...)
  • Save MODFLOW grid with the appended data to a Pickle

Plotting function to check that raster data was correctly applied to the MODFLOW data. Uses a colormap so colors cannot be manually defined

Inputs:

  • a MODFLOW geodataframe that was appended with data i.e, a 'value' column

Outputs:

  • a cloropleth plot of the MODFLOW grid with the appended data

Improvements:

  • Currently uses LUid to mask out empty/nan pixels. This could be fixed so that it works for any raster
  • Change LUid for something more general (like 'value')
Loading output library...

Plotting function to check that raster data was correctly reclassified and mapped to the MODFLOW grid. Uses a colormap so colors cannot be manually defined

Inputs:

  • the original raster
  • the reclassified raster
  • the MODFLOW grid appended with the raster values

Outputs:

  • 3 suplots comparing the datasets

Improvements:

  • Currently uses LUid to mask out empty/nan pixels. This needs to be fixed so it works for any raster
  • Change LUid for something more general (like 'value')
Loading output library...
Loading output library...
Loading output library...
Loading output library...