Thursday, August 13, 2015

Spinning-up a watershed model

by Laura Condon and Lauren Foster
Keywords: Domain setup, watershed model, getting started

Manual Location
Spinup options: 6.1.34
Directory of test cases: 3.5
Richards Equation Solver Parameters: 6.1.33
Boundary Conditions Pressure: 6.1.24
Manipulating Data (PFTools): 4

What is spin-up and why do it?
After assembling your ParFlow domain you will want to initialize your model to achieve some form of dynamic equilibrium before doing other testing on your domain. This is an important step for two reasons:

(1)         If you start running tests before your groundwater has equilibrated, then the changes you are looking for in your simulation can be obscured by the evolution of the system.
(2)         Initializing the model is a good way to get your domain solving smoothly before you run your simulations.

Of course there is not ‘right way’ to do this and its up to you to determine the best approach and what metrics you will use to determine that your model is adequately spun-up depending on what you plan to use it for. Refer to the reference list below for some articles on spin-up and the importance of initial conditions.

A general approach that we have used is as follows:
  • Initialize your water table. You can experiment with this but common approaches are to either start with a constant water table depth, start with the domain completely dry or initialize the water table to some realistic guess of how you think the groundwater should be configured.
  • Run for a long time (probably hundreds of years) with a constant recharge forcing or repeating some historical transient forcing until your water tables and groundwater storage stop changing. Usually you can accomplish this step with the surface flow turned off at first to make things easier to solve (see following sections for how to do this). 
  • After the groundwater is at some steady state or quasi steady state you can turn the surface flow back on and run until all of your streams have formed and again you don’t see changes in your groundwater storage. If you are planning on running ParFlow-CLM it is a good idea to spin-up the model without CLM first. This way you can get the groundwater surface water systems running smoothly before you turn on all of the land surface processes. Running this simplified version will speed up the process.
  • If you are running with CLM the next step is to turn on CLM and run for several years until the changes between years are small. Note that these are transient simulations so you are never going to achieve a steady state solution. You just want to be sure before you start running test cases that you are not gaining or loosing significant volumes of water to the subsurface over time (unless of course this is what you are expecting to see).

More specific instructions for setting up spin-up runs are provided below. Also, the LW_var_dz_spinup.tcl test case (you can find this in the test directory where you installed ParFlow) is an example where the water table is initialized at a constant depth with a constant recharge forcing is applied. 

Initializing the water table:
You can set the water table to a constant depth like this:
pfset ICPressure.Type                                   HydroStaticPatch
pfset ICPressure.GeomNames                              domain
pfset Geom.domain.ICPressure.Value                      -10.0
pfset Geom.domain.ICPressure.RefPatch                   z-upper

or start your domain off completely dry like this:
pfset ICPressure.Type                                   HydroStaticPatch
pfset ICPressure.GeomNames                              domain
pfset Geom.domain.ICPressure.Value                      0.0
pfset Geom.domain.ICPressure.RefPatch                   z-lower

or you can read in an input file with pressure heads like this:
pfset ICPressure.Type                                   PFBFile
pfset ICPressure.GeomNames                        domain
pfset Geom.domain.ICPressure.FileName

Setting up a constant recharge forcing:
If you are running without CLM (usually advisable until you have a stable water table and streams have formed) you can set a constant flux over the entire top of the domain like this (see manual 6.1.24):

pfset Cycle.Names                                  "constant"
pfset Cycle.constant.Names                    "alltime"
pfset Cycle.constant.alltime.Length                   10000000
pfset Cycle.constant.Repeat                              -1

pfset Patch.z-upper.BCPressure.Type                  OverlandFlow
           pfset Patch.z-upper.BCPressure.Cycle                  "constant"
pfset Patch.z-upper.BCPressure.alltime.Value       -0.0001

Note that the flux value is negative because you want water to flow down into your domain, which is the negative z direction. Also note that the units are [L/T].

If you want a constant flux in time but you want it to vary spatially you can do this by creating a pfb file that has the fluxes for every grid cell in your domain like this (see manual 6.1.33)

pfset Solver.EvapTransFile            True
pfset Solver.EvapTrans.FileName  "PmE.flux.pfb"

This input file should be 3D so you will need to fill in the vertical layers that you aren’t applying the flux to with zeros. Also, the units of this are not [L/T] they are [1/T]. So if you have a recharge forcing that is in [L/T] you need to divide by the dz of the layer you are applying the flux to get into the correct units.

Turning off the lateral flow during spin-up:
A common approach to spin-up is to start by turning off lateral overland flow so that ParFlow doesn’t have to solve the coupled surface-subsurface system while the groundwater is reaching a dynamic equilibrium.  You can do that as follows:
pfset OverlandFlowSpinUp           1

This key removes any ponded surface water at every time step so that no overland flow occurs. This is just a trick to make the problem easier to solve in the beginning and should not be used for regular simulation. After it appears that the water table is no longer changing you will want to set this back to zero and keep running to let all of your rivers form normally.

Other keys to speed up spin-up:
There are also two dampening keys detailed in section 6.1.34 of the manual. These keys make the problem easier by dampening pressure relationship in the overland flow equation with the following dampening term that gets added to the equation when pressures are less than zero:


 Values of P2 and P1 can be set as follows:

pfset OverlandSpinupDampP1  10.0
pfset OverlandSpinupDampP2  0.1

Note that the values of 10 and 1 for P1 and P2 respectively are just examples. You can adjust these values as you like to improve your solution.

Again these keys are useful for getting the problem started, but once you have things behaving nicely you should turn these off and they should not be used for regular simulation.

In addition to the spinup flags you may also need to make your timesteps smaller or adjust your solver settings. This is really a process of trial and error and it will vary depending on what your domain looks like. For tips on monitoring the solver performance and trouble shooting slow models refer to this blog post.

How do you know when you are done?
As noted above how you approach spin-up and how you decide when you are done is a subjective matter and it will vary depending on the questions you want to answer with your model. As you are spinning up some good things to look at are: changes in total subsurface storage, total surface storage, total overland flow exiting the domain and the flow at a point of interest. You can calculate these using your own scripts based on the and out.satur files or you can use the built in PFTools. Examples for how to calculate each of these items are included in the manual section 4.3 (examples 6-9).

If you are spinning up with a constant recharge forcing you can check that there is no longer any change in groundwater storage either by calculating the total groundwater storage at two points in the spin-up or by subtracting your total model outflows from the inflow.

If you are running transient simulations with CLM the groundwater storage at the beginning of a water year should be very close to the storage at the end of the water year. This means you need to run your model for at least one full year. Calculate storage at these two time steps and compare. It may also be more useful to compare your change in storage to total precipitation input in your domain. At the end of the day the most important thing is that the order of magnitude of the result you want to look at is significantly greater than the change in storage over a year, otherwise it’s impossible to tease out your results from the uncertainty in the model spin-up.

Ajami, H., McCabe, M.F., Evans, J.P. and Stisen, S. Assessing the impact of model spin-up on surface water-groundwater interactions using an integrated hydrologic model. Water Resources Research, 50, doi:10.1002/2013WR014258, 2014.

Ajami, H., Evans, J.P., McCabe, M.F.and Stisen, S. Technical Note: Reducing the spin-up time of integrated surface water–groundwater models. Hydrol. Earth Syst. Sci., 18, doi:10.5194/hess-18-5169-2014, 2014

Seck, A. Welty, C., and Maxwell, R.M. Spin-up behavior and effects of initial conditions for an integrated hydrologic model.  Water Resources Research, 51, doi:10.1002/2014WR016371, 2015.

1 comment:

Surajit Deb Barma said...

I am wondering if spin-up studies on a large scale river basin would be possible for Ganges-Brahmaputra river basins (subsurface, surface/land surface interaction) with an area of 1.7 million km2. I would like to carry out my PhD project using ParFlow-CLM using ancilliary and observed data sets.

Our institute is in the process of acquiring HPC/ super computer and it might take a while to have one. Hence, I am planning to use national facility center via remote access. However, many researchers might use at the same time. Could you suggest me how long it might to run 1 year of simulation at 1 km horizontal resolution (vertical resolution not yet decided) by using atleast 20 % capacity of the super computer (details attached)? I have mentioned the details here as well.

PARAM YUVA-II national facility has a peak performance of 529.74 TF (Teraflops) and Sustained performance (Linpack) of 383.708 TF with a total nodes of 225. Do you think it would be feasible to carry out such studies given huge computational requirement?

Thanking you.


Surajit Deb Barma