Running pywatershed with time_step > 1 day

[cneyens]

I'm trying to run pywatershed on a weekly time step. I've resampled the meteorological forcings to weekly basis (summing for precipitation, averaging tmax and tmin) and I've set the time_step argument to a time delta of '7D' (and changed the time_step_unit option to 'D'). It seems pywatershed can handle this, but the swrad values are not adjusted properly, resulting in an underestimation of the evaporation fluxes compared to the run with a daily time step. I've been trying to use rad_adj to multiply swrad by the time step length, but to no avail. Also, far less interception is occuring in the weekly model. I guess some of the parameters need adjusment, but it is not clear to me which ones.

Is there a way to have pywatershed run with a time_step > 1 day? This would reduce the runtime of our large 30-year model substantially, since we're only interested in groundwater recharge on a weekly to monthly basis. If not, what is the role of the time_step parameter in the control file?

[jmccreight]

Currently, the timestep is just an ambition. There are a few other places this shows up, but none are actually functional. I'm sorry this is misleading. There is really only daily timestepping available (with some channel flows being calculated hourly). Moving to an arbitrary timestep is not really a priority that has come up.

Yes, and the parameters would all need recalibrated, as you mention.

[jmccreight]

Speaking of runtime, I'm working on improving run times. This would likely include the ability to reduce output resolution as well, in time. This is definitely many months out, however.

Just to be sure, you are using calc_method="numba" where available? Are you also turning off processes you dont need. As far as I can infer, this would only be the PRMSChannel in your case of wanting recharge.

[cneyens]

Thanks for the reply!
Yes, I'm using numba and I'm turning off PRMSChannel, PRMSSnow (negligble in the area of interest) and PRMSGroundwater. The long runtimes are mainly due to the gridded HRU approach we're using (50x50 m).

[jmccreight]

Oh, yes, I remember now that you'd turned off Snow and are using gridded.

What is the total number of spatial elements/HRUs? You might get some speed up by configuring multithreaded numba options. I've seen this in some gridded cases with number of elements > ~8k. Here's the code to put above your imports

import os

os.environ["NUMBA_NUM_THREADS"] = "2"
os.environ["NUMBA_THREADING_LAYER_PRIORITY"] = [
    "omp tbb workqueue",
    "tbb omp workqueue",
    "workqueue omp tbb",
][0]
calc_method = "numba"

I think the "NUMBA_THREADING_LAYER_PRIORITY" is pretty platform dependent. You might try all of them on some shorter period of time. I've never found going about 2 threads to improve much, but maybe if your domain is really big it might help. The NUMBA_NUM_THREADS env var is used to set the number of threads for each process, nb_parallel. From this search, you can see that it appears that Canopy, Snow, Runoff, and Soilzone are using it (and Groundwater is not). LMK if you have any success with that.

I think the bigger problem is that the time loop leaves numba and re-enters the python interpreter between processes. This is what I'm trying to reduce/minimize in the bit of work I mentioned previously: keeping periods of time of calculation out of the interpreter.

You've probably also thought of this, but what is the correlation length scale of your recharge? you could drop to that resolution and then do a simple downscaling/interpolation to 50m.

[cneyens]

We are basically creating a tool to estimate groundwater recharge for use in local groundwater modelling studies by wrapping pywatershed + MODFLOW6-UZF (one-way coupling at the moment, so no API). The end-user can define an area of interest, a time-period interval, and the output resolution. Base input is provided at 50x50 m resolution for the entire region (climate forcings at 5x5 km), but the user can provide their own lu/soiltype maps + parameter sets or forcings (which would make the downscaling tricky since the correlation lengths may vary depending on the input).

I'm running a test-case of 10x10 km for 6 years, so about 40k HRU's. Right now, everything works as intended and the results are there, with runtimes of about 7 min on my machine. For a larger area/longer period, I expect these to become excessive. I've tried the numba options but they didn't seem to be doing much in terms of runtime. We were thinking of using Dask to split up the domain into chunks and running those pywatershed models in parallel on a cluster. Since cascading flow is inactive, all HRU's are effectively independent of each other, but perhaps there is little gain in runtime compared to the numba approach.