SalvIac / pysimulator

Simulator package for earthquake catalogs - Salvatore Iacoletti
GNU Affero General Public License v3.0
2 stars 1 forks source link

Incorrect time-dependent incompleteness #2

Open mherrmann3 opened 2 years ago

mherrmann3 commented 2 years ago

I found that your implementation of the incompleteness model still leaves some events below the completeness threshold. I found the possible reason:

Instead of (jump to line)

np.clip(compl_vs_time_general(mag_par, deltat, *coeffs), mmin, mag_par)

I'd propose:

np.maximum(compl_vs_time_general(mag_par, deltat, *coeffs), mmin)

But this does not treat cases with overlapping aftershock sequences (e.g. when a large event like the mainshock dominates all other "incompletenesses"). For this reason, I think it's better to do a final "cleanup" by looping again over all events and applying this constraint, like so using an "Mc mask":

mc = np.ones(len(stcat_df)) * params['min_mag']

for i, row in stcat_df.iterrows():

    rel_time = (stcat_df.datetime - row.datetime)  # unit: days
    mask = rel_time > 0
    mc[mask] = np.maximum(compl_vs_time_general(row.magnitude, rel_time[mask], *coeffs), mc[mask])

df_withSTAI = pd.concat([stcat_df, pd.Series(mc, name='Mc')], axis=1).loc[stcat_df.magnitude > mc]

(that's how I applied STAI to the returned catalogs).

Thanks again for all your hard work, Sal!

SalvIac commented 2 years ago

Nice! I'll look into this.

Do you mind sending a small example where you've noticed that np.clip doesn't work?

Also, on the overlapping sequences: I think the issue here is more than just "coding". Available magnitude incompleteness models don't have a spatial component (a larger mainshock might dominate a specific area while another dominates an adjacent area). My version of the clipping assumes that every aftershock is influenced by its mainshock only. This is surely an approximation, but I think the spatial component is also needed to apply more sophisticated models with different mainshocks affecting an event. Of course, this only kicks in if two very large mainshocks are generated in a super short time (i.e., less than a day usually), super close to each other. So I am not sure it drastically changes the results, but it will be worth including if/when someone comes up with a model for a time-space dependent magnitude incompleteness model.

Thanks again for the helpful suggestion. Sal

mherrmann3 commented 2 years ago

Of course I can provide an example! With the parameters of example 1 and "incompletess:True", I simulated a catalog for a Ridgecrest-like sequence:

example1_simulation+incomple-withSTAI

The time-varying completeness approach I proposed in the initial post ("Mc mask") is calculated at each event and shown in red; black dots are synthetic events above this threshold; gray ones below.

Because I didn't need the internal incompleteness functionality anymore, I set "incompletess:False". Then, pysimulator outputs a few more events (324 vs. 265), but with a completely different magnitude behavior:

example1_simulation-withSTAI

It appears as if "incompletess:True" does not only cut events away, but shifts them above the time-dependent completeness level (still keeping some below), creating much stronger aftershocks.

I'm actually not sure whether "incompletess:True" or "incompletess:False" behaves wrong (I guess the former, because those large aftershocks should create much more productive secondary aftershock sequences), but I would expect that the distribution of magnitudes above the red curve is about the same in both cases (visually).