Prevent numerical instability in pt_downscale_radiations

[mmazzolini]

Hi everyone,

I have experience some numerical instability when downscaling the Shortwave Radiation. I am downscaling ERA5 to a set of representative points in a basin. Some of them have large slope and I think that's why I encountered this problem.
I noticed the problem as in some points I got nonphysical values of SW, caused by a large SW_direct component. Here's an example of a west-facing and 30 degrees steep point:

I get a spike of SW radiation in the evening (sorry for the wrong time axis, times are in UTC but the basin is at UTC-7). That is caused by the fact that the ratio between cos_illumination and mu0 explodes, likely because there is a very small value of mu0 just before the sunset.

I thought of solving it by setting a minimum value of mu0 as 0.25 mu0_stable = mu0.where(mu0 > 0.25, 0.25). That value corresponds to a sun zenith of 75 degrees. I called this mu0_stable and substituted it in the formulas where mu0 appears in the denominator.

Moreover, while I was exploring these terms I also discovered that, for some approximation errors I guess, the calculation of surf_interp['SW_direct'] leads in some cases to values slightly larger than the SW radiation on the top of the atmosphere (SWtoa). When that happens, I got that ka was at times negative, as you can see in the next figure (representing another west-facing point but steeper). Negative ka also leads to exploding SW_direct, this time in the formula for down_pt['SW_direct_tmp'].

So I removed the values of surf_interp['SW_direct'] that are larger than SWtoa by substituting them with SWtoa:

surf_interp['SW_direct'][surf_interp['SW_direct'] > SWtoa] = SWtoa[surf_interp['SW_direct'] > SWtoa]

I've already discussed about the first issue with @krisaalstad. However, I'm not sure at all that these fixes work well in other cases. I guess there needs to be some testing in other basins before accepting these?
If needed I can provide you with the dataset I'm working with.

Cheers,
Marco.

[ArcticSnow]

Thank you @mmazzolini ! Good catch.

I won't have time soon to verify this non physical spikes occur in for instance the the example watersheds. Or in case you have time on your hand, I'll let you check too with the example cases: https://github.com/ArcticSnow/TopoPyScale_examples/tree/main

Also, we now have another downscaling workflow where this would need to be corrected too: https://github.com/ArcticSnow/TopoPyScale/blob/main/TopoPyScale/topo_scale_zarr.py I put this for when we follow through the PR to not forget to correct the bug here as well, if consistent.

This points to the need to put in place some thorough testing scripts, looking at hourly behavior of the output variables.

[joelfiddes]

Hi @mmazzolini , thanks for this catch and sorry for slooow response!

The fix looks reasonable overall. A few suggestions (for all of us) before merging:

1. Hard coded constant
The 0.25 threshold for mu0 should be defined as a named constant with a brief comment explaining the rationale (sun zenith ~75°, beyond which the corrections become unreliable).

2. Use .clip() or .where() instead of in-place assignment
This pattern can have unexpected side effects with xarray:

surf_interp['SW_direct'][surf_interp['SW_direct'] > SWtoa] = SWtoa[...]

Cleaner would be:

surf_interp['SW_direct'] = surf_interp['SW_direct'].clip(max=SWtoa)

or equivalent with .where().

3. Root cause of SW_direct > SWtoa?
The cap is pragmatic and physically correct, but do you have a sense of why this occurs? Interpolation artefact? Elevation correction? Might be worth a follow-up issue to investigate even if we merge this fix now.

4. Zarr workflow
Need to include the equivalent fix in topo_scale_zarr.py in this PR? Otherwise we'll have divergent behaviour between the two pathways.

5. Test case
Would be great to add a minimal test case - a synthetic steep west-facing point that would have triggered the bug pre-fix. Doesn't need to be elaborate, just a few timesteps.

Thanks again for tracking this down :)

[joelfiddes]

@ArcticSnow @mmazzolini in general I think we are good to merge subject to comments above as both are “no regret" solutions.

thanks again!

[joelfiddes]

@mmazzolini @ArcticSnow just had chat with Simon - I can handle the PR from here if you like, would still be great to have link to a minimal example for reference, particularly in eventually tracking down the root cause. cheers all