Model spinup refers to the pre-production run period that you run your model for to let things equilibrate before you start testing your hypotheses. In this post, I also include some strategies I used to get my model up and running.
I am applying ParFlow to a very small domain (0.05 km2) to model fine-scale urban hydrological infiltration and runoff processes. My horizontal resolution is 1.5 m x 1.5 m, one of the finest scale grid resolutions that I have seen ParFlow applied to (others in my my lab are also working at these fine resolutions). At this scale, the surface and subsurface dynamics of overland flow and infiltration/movement through the subsurface can happen at vastly different order of magnitude, often resulting in non-convergence of the solvers at the core of ParFlow.
I went through three distinct spinup “phases” that eased some of the non-convergence issues I was seeing.
- Start with homogeneous permeability in the subsurface (run for 100 days at 1 hour timestep)
- Add in permeability heterogeneity of subsurface and add constant recharge forcing (run for three years at 1 hour timestep)
- Add in NLDAS 1D meteorological forcing (run for 6 years at 1 hour timestep)
In Step 1, I started with an initial pressure condition with the water table set 4 meters from the surface (knowing that historical depth to water in my domain area ranges from 4.88 m to 6.25 m). I set the four side “wall” boundary conditions at constant pressures of 5 m from the top face (I am using a Terrain Following Grid setup). This means that equilibration over time will be draining the domain from all sides so that the water table is falling towards 5 m from the surface, and the falling water table will be propagating from the four side walls of the domain. I am not interested in overland flow, so I turn on both spinup and dampening keys in my tcl script.
In Step 2, I distribute my 3D permeability pfb file that includes the specified permeabilities for various zones within my domain (paved surfaces, buildings, pervious surfaces, etc). I also add in a constant hourly boundary condition to the top of the domain. This is the annual cumulative recharge to the subsurface of about 1′ that I got from a geologic report of the region. No other processes (evapotranspiration, etc) are included at this point and spinup and dampening keys are still on. After three years, there is still about 0.25% change in subsurface storage in my domain.
Change in subsurface storage for Step 2:
Year 1 to Year 2: -3.5%
Year 2 to Year 3: -0.68%
Year 2 to Year 3: -0.28%
In Step 3, I replace the constant forcing with NLDAS meteorological forcing data, which includes 8 parameters (precipitation, surface pressure, temperature, 2 directions of windspeed, humidity, and 2 directions of radiation). These parameters are used in the CLM portion of the model, which also takes as an input the vegetative surface cover to calculate spatially varied evapotranspiration processes.
Change in subsurface storage for Step 3:
Year 0 to Year 1: -0.300%
Year 1 to Year 2: -0.236%
Year 2 to Year 3: -0.168%
Year 3 to Year 4: -0.126%
Year 4 to Year 5: -0.104%
Year 5 to Year 6: -0.087%