Augmenting rupture rates for classical PSHA with aftershocks

[cossatot]

We are working on an implementation of classical PSHA with aftershocks.

The current approach that we are taking is to parse the source model, and place all of the aftershocks for each mainshock rupture on the other ruptures already in the SSM, yielding an additional rate increment for each rupture, conditional on the occurrence of the mainshock ( see https://gitlab.openquake.org/hazard/projects/2020/metis/-/blob/master/t442/openquake/aft/aftershock_probabilities.py). These rate increments are then converted to unconditional rates, and then summed for all of the mainshock ruptures. Finally we have an extra rate increment that describes the rate at which each rupture in the SSM will act as an aftershock. This value is to be added to the primary (mainshock) rupture rate before the calculation of the PoEs during a typical classical PSHA. All of the steps before the addition of the aftershock and mainshock rupture rates are performed during a pre-processing stage, with code that for now won't be included in the Engine (among other reasons, it introduces Numba as a dependency; whether to include it in the future is a separate discussion).

The additional rupture rates are stored in a CSV file, with a format like:

rate,source,rup_num
0.00045937636367182145,i17,0
0.00046102771506415717,i17,1
0.0004443637497429468,i17,2
0.00046275940241261834,i17,3
0.0004931585637198938,i17,4

where source is the source_id of the source that the rupture belongs to, and rup_num is the number of the rupture when the source is enumerated with iter_ruptures().

There does not need to be any configuration parameters for the job.ini except for the file path to the CSV file. The necessary scientific parameters are used during the pre-processing step. So the job.ini just needs an entry like:

[aftershocks]
rupture_occurrence_rate_adjustment_file = ./rup_adjustments.csv 

I have a very crude implementation in the Engine for reference here: https://github.com/gem/oq-engine/compare/rup_prob_adjustment. The code reads the rupture rate adjustments file, and then when any source is processed in openquake.hazardlib.contexts.ContextMaker.gen_ctxs(), the occurrence rate for that rupture is adjusted just before the contexts are created (https://github.com/gem/oq-engine/blob/99ae182965ddeb7d22e978eb365c9507650a0d3a/openquake/hazardlib/contexts.py#L403). This was the way that made the most sense to me, and it works fine except that because I don't understand how the sources are partitioned when they are split for parallel processing with Starmap, I had to turn off the source splitting (so that the ruptures can be correctly indexed by the rup_num), so it is quite inefficient. I don't expect that this code will be used in production, it is just a crude reference implementation.

I am happy to answer any questions and provide source models and rupture adjustment files for testing.

Thanks!