STEPS-blending: Losing small scale features in the first time steps and other issues

[mpvginde]

Hi everyone,

In setting up pySTEPS for a blended nowcast with STEPS, I have noticed that during the first time steps some of the radar image small-scale features are immediately lost.
This becomes particularly apparent when the blended STEPS nowcast is compared to the non-blended version of the STEPS nowcast.

An example over the Belgian radar domain can be seen here:

Fig. 1: Identical starting conditions (i.e. radar field) for the blenden nowcast (left) and non-blended nowcast (right)

Fig. 2: The blended (left) and non-blended (right) nowcast after 1 time step of 5 minutes

It's clear that even after the first time step (Fig, 2), the blended nowcast field has lost some fine-scale structures w.r.t. the non-blended nowcast (i.e. the blended fields seem smoother). In addition, there also seem to be more broad structures with more intense rainfall in the blended nowcast.

When investigating this phenomenon I came across several issues/points of attention in the blending code.
I already discussed some of them with @RubenImhoff, but I will try to explain them here too so everything is well documented.


In order to explore the problem I plotted the time evolution of the Power (displayed as standard deviation [St. Dev.]) of the different cascade levels for the different fields.
For blending nowcast there are 2 relevant fields per level:

• the precipitation field
• the noise field

Both are evolved separately using an AR(2) model
In addition I also plotted the time evolution of their respective blending weights with which they are combined. I have deliberately omitted the NWP-component as its weight is typically very small during the first time steps.

For the non-blending nowcast there is only 1 relevant field as the noise is directly incorporated in the AR(2) process.

Finally, also included is a figure of the actual fields themselves at a specific time step.
For the blended nowcast:

• the AR(2)-evolved precipitation field
• the AR(2)-evolved noise field
• both fields combined using their respective blending weights

For the non-blended nowcast:

• the AR(2) with noise term-evolved precipitation field

Fig. 3: AR(2)-evolution overview for cascade level 6 at time step 1

Fig 3. shows this AR(2)-evolution overview for cascade level 6 (2nd smallest level).
We notice some problems:

1. The AR(2) evolution of the precipitation cascade without the noise term EPS is no longer stationary in the variance (something that is necessary according to Bowler et al. 2006:
   "The noise term has been omitted from this formulation and the cascade is renormalized at each step to ensure no loss of power")
   In practice, this means a dampening of the amplitudes (= loss of power) for all cascade levels
2. The noise starts from 0 and needs to spin-up. For the largest levels (due to the large temporal autocorrelations) this may take a long time.

This results in the following issues:

• The precipitation component gets punished twice: Once because the power is decaying due to the AR(2)-evolution and once due to the weighing by the blending weight
• Also the noise amplitude is not high enough because the noise needs to spin up (started from zero) and is also weighted by the blending weight.

What does this mean:
For the smallest levels, due to the low temporal autocorrelation, the precipitation field decays very fast and also the blending weight immediately is very low. This mean that the original precipitation signal is almost immediately lost. This can be seen by comparing the blended field from the blended nowcast with the field from the non-blended nowcast (the two right-most panels in the bottom row of Fig. 3).
Additionally also the noise has a reduced amplitude (due to its initialization to zero) which results in a blended field with a too small power (= St. Dev.) as seen in the right-most upper panels in Fig. 3.

For the largest scales, there is also a loss of power in the blended field. This is mainly due to the spin-up of the noise field. Due to the large temporal-autocorrelation they stay close to zero for a long time. Keeping in mind the blending strategy, I guess the St. Dev. of the noise field should be (close to) 1 and the blending weights should handle their contribution to the blended field. But since this is not the case the overall power of the blended fields is too low for the largest scales (as seen by comparing the two right-most upper row panels in Fig. 4).

Fig. 4: AR(2)-evolution overview for cascade level 2 at time step 10

These issues lead to two effects: The first one being the loss of small scale features in the first time steps of the nowcast (which triggered this investigation).
The second one, which is more subtle, being the large scale noise having not enough power in the beginning of the nowcast.
This syndrome can be seen when there is precipitation from the NWP-forecast outside the radar-domain (e.g. in Fig. 2). In this case the noisy fields (since the NWP blending weight is typically rather small) appear with rather small rainfall amounts (which actually gives them a quite realistic feel, but theoretically their amplitudes are too small.)

Unfortunately there is another problem. Namely with the semi-Lagrangian advection of the small scales.
In the blending nowcast all the different cascade levels need to be advected separately because they need to be blended with the NWP cascade levels. But for the smallest levels this also leads to a loss in power (as seen in Figs. 5 and 6)

Fig. 5: Standard deviation (~ power) of the noise field before and after semi-Lagrangian advection for level 6.

Fig. 6: Standard deviation (~ power) of the noise field before and after semi-Lagrangian advection for level 7.

The reason for this is less clear to me, but I suspect this has to do with the interpolation needed to find the starting points of the advection. Due to the high spatial variance of the smallest levels, potentially a large error is introduced when interpolating.
The non-blending nowcast does not suffer from this because their the levels are recomposed to the total precipitation field before the advection takes place.

Some possible solution:

1. Renormalize the precipitation cascade after each AR(2) iteration
2. Initialize the noise cascade with the appropriate noise

I did a first implementation if these ideas which resulted in the following:

Fig. 7: AR(2)-evolution overview for cascade level 6 at time step 1 after adaptation

Fig. 8: AR(2)-evolution overview for cascade level 2 at time step 10 after adaptation

As Figs. 7 and 8 show, the blended cascade levels now have again the correct St. Dev. (or power) after the AR(2) evolution.
The resulting nowcast is shown in Fig. 9

Fig. 9: The blended (left) and non-blended (middle) and adapted-blended (right) nowcast after 1 time step of 5 minutes

As you can see the adapted nowcast retains better the small-scale features in the first time steps. Also notice the precipitation outside the radar domain immediately coming in at full power (a result of the noise component of the large scales having the correct power).

However, these adaptations still have some drawbacks:

• renormalizing the precipitation field after the AR(2) evolution results in a constant precipitation field. So this technique is not ideal. An option could be to use the localized version of the AR(2) process, but this becomes quickly very computationally intensive
• No fix is found for the loss of power in the small scales due to the semi-Lagrangian advection. This is also seen in the power spectrum shown in Fig. 10.
  
  Fig. 10: Power spectrum of the radar field, blended nowcast, non-blended nowcast and adapted blended nowcast after 1 and 2 time steps. From this graph it is clear that even after the first time step the blending nowcast loses power at the small scales, while this is not the case for the non-blending nowcast. The adapted blending does perform better but there is still some loss of power in the smaller scales.

@ladc and myself have reserved some time in July to further investigate this issue.
Comments and suggestions are much appreciated.

Kind regards,
Michiel

[RubenImhoff]

Hi @mpvginde,

Many thanks for this thorough analysis! We already discussed it earlier, but I'm very happy to see you went through the code and these results in so much detail. And good to see this as a pysteps issue now, so that we can think along in coming up with a solution.

Regarding initializing the noise cascade with the appropriate noise, I think this is the way to go and what you have tested, could in my opinion directly be implemented in the code.

Regarding the renormalization, I remember that @ladc and I discussed this and I can't entirely remember why we decided not to implement that (directly). I think we were worried about any issues at the end of the forecast when the x times renormalized cascades had to be recomposed. Have you tested this by any chance? If there are no issues, I think this would be a great addition as well. The localized version of the AR(2) process would even be better, perhaps this is eventually the way to go for the blending procedure. Maybe we should consider this as step 2 in the improvements? I can imagine that the localized version also takes somewhat more implementation time.

Finally, regarding the semi-Lagrangian advection related issues, perhaps @pulkkins has ideas here?

Cheers,

Ruben

[dnerini]

Very interesting issue and a lot of detail to digest! Many thanks for the analysis @mpvginde!

The noise starts from 0 and needs to spin-up. For the largest levels (due to the large temporal autocorrelations) this may take a long time.

I agree this doesn't sound right. Should we consider it as a bug or does this follow the original formulation in the paper? Can you submit a PR with the fix you've tested above?

The AR(2) evolution of the precipitation cascade without the noise term EPS is no longer stationary in the variance (something that is necessary according to Bowler et al. 2006

Again, should we consider this a bug? From you answer @RubenImhoff , it seemed that it was rather a design choice, was it?

Finally, I'd like to ping @albertocarpentieri here, since he encountered a similar issue in his work (preprint ), although it was related to the extrapolation only case without blending. But he also added a renormalization step after each AR iteration, so maybe there is a more general solution that we can discuss with respect to this?

[dnerini]

@mpvginde I sent you an invitation to join the pysteps organization. Feel free to create new branches directly in the pysteps repo.

[RubenImhoff]

Hi @mpvginde, is it also getting time to merge our blending improvements to the main branch of pysteps, as that will solve this issue?