We need some GIS data sets to get us started. We'll need a DEM, flow lines, and a watershed boundary to start. We'll need more as time goes on, but these are needed to get us started. I've provided started data from USGS. All of this data is in NAD83, geographic coordinates. It needs to be project to a flat earth. We can do this reprojection a variety of ways, but let's play with using python to do it. Check out ReprojectWithPython.pynb.
This script resamples the DEM to a desired grid for simulation.
- After we run this script, we'll have a resampled dem. However, the resampling is going to create artifacts and needs to be processed for surface flow calculations. We'll do this using the SAGA tool kit in QGIS.
- Processing Steps in QGIS using SAGA:
- Fill Sinks -> needs: resampled dem; produces: filled dem
- Use the planchon/darbou - minimun slope set to 0.01
- Catchment Area -> needs: filled dem: produces: upslope accumulated area
- Use recursive, D8, no manual sinks.
- Channel Network -> needs: filled dem, upslope accumulated area (for initiation points); produces stream network (raster and shape)
- Experiment with initiation threshold. 5e5 is a reasonable minimum starting point I've found.
- Burn Stream Network
- epsilion 0.5, method [1] -> needs filled dem, stream network raster, flow direction raster
- Fill Sinks -> needs: resampled dem; produces: filled dem
- Does some further reprocessing of the DEM and then plots flow accumulation and flow directions using the GSFLOW utilities, to get everything ready for modflow/gsflow.
In order to get watershed boundaries for our new raster, we'll need to find appropriate pour points. This is best done using QGIS. Open up a map fill your filled, burned DEM, along with the Rattlesnake Flowlines shape file and the original Rattlesnake Watershed shapefile. So, for some context, first of all the NHD flowlines that make up the Rattlesnake flowlines aren't perfect. They were our best guess at the channel sometime in the past. Riley knows that channels move. They probably weren't correct everywhere to begin with. But now, our DEM is not really the same as the real world any way, so we're going to have to do some sleuthing to find the correct pour points.
- Create a new layer: Layer->Create Layer->New Shapefile Layer
- Call it Rattlesnake.shp or something
- Make the geometry
point
- Make its coordinate system UTM 12N
- Save it in your
data/gis/
folder - I add a field called
Drainage
where I will store an identifying name.
- Find the pour point.
- use the information pointer tool, to identify the lowest elevation cell at the mouth of the drainage.
- hint: I change the symbology of the DEM to pseudo color, and rescale the min-max to the current canvas after I zoomed into the area near the mouth.
- Edit the shapefile
- Add a point in the middle of cell of lowest elevation.
- label the drainage "Main Rattlesnake Creek" or something like that
- Save edits
- Create another shapefile called RattlesnakeSubSheds.shp
- Using the same procedure as above, add pour points for all branching drainages seen in flow lines (except the irrigation ditches).
- The pour points need to be above the main Rattlesnake channel, but as close as you can get to it. It takes some clicking around to figure out what's going on.
- You can figure out the name of the subwatershed for labeling by using the information pointer on the flowlines shapes file.
- Hint: you might want to rescale the symbology occasionally to help visualize elevation differences near a given confluence.
This script calculates the uplope areas for all your different pour points. It it needs the name of your Rattlesnake pour point shape file made above and the subsheds shapefile name made above. It should produce a figure with the different watershed areas, as well write some file for our use later.
In this script, we'll create stream network and cascades, which are used in a subsequent step to define the location of stream cells in the model for the stream boundary condition.
- Make sure to the file path names to match your directory structure.
- Much like creating stream networks using SAGA, we'll need to play around with an area threshold for initiating the stream. This variable is called
threshold_m2
in the script. Try to match the NHD flowlines as best as possible.
In this script we'll create in initial MODFLOW model. This will be a super basic model, just to get a forward MODFLOW model running. This model will then be modified to create more complicated versions. To get this model running, you mostly need to:
- Get paths all managed up
- In particular, the executable needs to be in the path.