Time-dependent super-droplet sampling and modification

[jtbuch]

Hi @slayoo,

I have been running experiments with PySDM to simulate the injection of a new 'seed' aerosol distribution with distinct properties than the background aerosol distribution used to initialize the simulation. As you might anticipate, this functionality has several applications under the cloud seeding research umbrella, including marine cloud brightening, in-situ growth of ice crystals, and precipitation enhancement. For the rest of the discussion in this issue, I will only focus on the precipitation enhancement case with two modifications of the default PySDM code:

1. Adding new super-droplets: In my first attempt, I tried simulating the injection of seeded aerosols at arbitrary time intervals by adding new super-droplets (SDs). This involved writing a new inject_particles (see here and here) method in the simulation file, making minor tweaks in relevant methods of the Displacement and Collision classes, while also writing a couple of small setter methods to keep track of the new SDs to avoid array shape mismatches.
   The relevant notebook for running the experiment in a SH2012 environment is available here, where I try to show that despite all the modifications, both displacement and collision-coalescence fail, i.e. there is no change in multiplicity and positions of the seeded SDs, even though the code executes successfully.

2. Modifying existing super-droplets through the 'hijack and redistribute' method: After several stimulating conversations with @claresinger, I implemented her suggested approximation of aerosol seeding: a) 'hijack' a cloud SD of a size similar to the desired seed and update the value of all relevant extensive attributes: multiplicity, kappa, water mass etc., b) 'redistribute' the original cloud SD's attributes to its neighboring SDs in a grid cell. This way, I can 'inject' a new seed SD with the appropriate multiplicity attribute to represent the integrated aerosol concentration.
   As the spectra, coalescence rate, and rain water mixing plots in the associated notebook indicate, there is a distinct seeding signature for well-chosen values of seed radius and concentration.

I'd be curious about your thoughts on this approach overall! Specifically, if there is a way to get the code working while adding new SDs, or if there is another way to model seeding in PySDM. One suggestion that both @claresinger and @edejong-caltech made (also something that @kdlamb and I discussed independently) was to introduce a new attribute vector for the SDs that acts like a flag to decide whether a particular SD participates in each step of the simulation. In such a setup, certain seed SDs will behave as 'ghosts' for a while before they are activated at the desired time step and then behave normally.

[claresinger]

Following up here after a discussion with @jtbuch last week...

My original suggestion was actually a bit different from what is described above as method 2. I'd call it something like “rebinning" the droplet distribution "via fake collisions”. The idea would be that for each timestep where you want to inject seed particles, you rebin the existing SDs in n-1 particles and use the newly freed SD to represent your seeds. The nuance is all in how you do the rebinning - you want to do it in a way that conserves the properties of the original distribution -- in obvious ways like mass and water conservation, but also preserving the character of the distribution, e.g. not getting rid of a long tail of few, large droplets.

Such an algorithm might look as follows...

For each time step when you want to inject seed particles:

1. Randomly sample from your seed particle property distribution to determine the attributes of the new SD you want to seed.
2. Search each grid box for the two most similar SDs. Include in this search your candidate seed SD (in case there is an existing SD that looks very similar to the seed). “Most similar” will be the tricky thing to define properly. But as a simple example it can be the euclidean distance in the n-dimensional attribute space (dry radius, wet radius, kappa). You should also include
3. Combine these two SDs in a “fake collision/coalescence” event, where the probability of coalescence is 1.
4. Update the attributes of the now “empty” SD with the properties of the seed particles, including multiplicity as determined by injection rate.

All the art here will be in defining the “similarity” metric. You want to ensure that your "rebinning" is conservative -- i.e. this procedure in the limit of very, very low injection rate does not change macro scale properties of the simulation (onset of precipitation, etc.). Maybe the distance in wet radius should be more heavily weighted than the distance in kappa, as one example.

[slayoo]

Thanks @jtbuch & @claresinger!

One suggestion that both @claresinger and @edejong-caltech made (also something that @kdlamb and I discussed independently) was to introduce a new attribute vector for the SDs that acts like a flag to decide whether a particular SD participates in each step of the simulation

In collisions and out-of-domain precipitation, we flag particles that should no longer take part in the simulation by setting multiplicity to zero or setting the index to n_sd:
https://github.com/open-atmos/PySDM/blob/main/PySDM/backends/impl_numba/methods/collisions_methods.py#L1117-L1128

This ensures that these are properly treated later on in condensation and coalescence summations. Perhaps there is no need for a new attribute?

[slayoo]

@jtbuch, @claresinger, apologies that it took longer than expected, but here's a first try at a new Seeding dynamic for PySDM based on the above idea of using zero multiplicities as flags for implementing a buffer of super droplets to be used for seeding: #1367

There is a draft of a simple parcel-model test case currently included in the PR: https://github.com/slayoo/PySDM/blob/sd_injection/examples/PySDM_examples/seeding/hello_world.ipynb

It includes condensation, coalescence and seeding. The seeding logic parameters are not yet robust to temporal resolution - there is just prescribed rate per timestep.

The output plot looks like this (blue lines correspond to background aerosol, red lines to the injected particles):

Feedback welcome!
In particular, it would be great to identify a simple test case from literature that we could refer to.

[slayoo]

(I've refactored a bit the PR today moving the plot-generating example into a Jupyeter notebook - hope this will facilitate trying it out: https://github.com/slayoo/PySDM/blob/sd_injection/examples/PySDM_examples/seeding/hello_world.ipynb)

[claresinger]

Hi @slayoo! Thank you - this looks like an awesome first step. I see that you have implemented a time window for the seeding, but not yet a variable for the injection rate, is that correct? And is it also correct that in this example the seeded particle attributes are fixed?

I think ideally we can have a specified injection rate and also inject particles from some specified distribution, but where the injected particles get randomly sampled from that distribution at each timestep. Maybe this means injecting multiple super particles, with different properties (dry sizes, kappas) and multiplicities at each timestep?

Possibly it will also interesting to be able to have more than one injection window: e.g. to be able to inject particles for 10 min every 1 hour... this could match up nicely with the variable updraft velocity from the Shipway&Hill examples. But I think if we can implement a variable injection rate then that can almost override the injection time window, by specifying the rate to be zero for some times and nonzero for other times.

@jtbuch How do these extensions sound to you? Are there other features you'd been thinking about for the experiments you want to run?

[slayoo]

@claresinger, thanks for following up. Yes, of course, the single time window, injection rate, single-super-droplet-per-step operation and lack of attribute sampling are all temporarily set here to the simplest possible settings. I was just aiming at clarifying how to technically address the idea.

[claresinger]

Yes, great! Just adding my two cents on what will be the interesting/physically relevant extensions to work on next. And hopefully my comments on the #1367 are helpful, too.

[jtbuch]

Hi @slayoo, I finally got around to running the example and reading the code in more detail. Your implementation of a buffer of super-droplets is quite elegant and much more robust than the strategy I have coded up.

Besides the comments raised above, I was wondering what the seeding would look like in space for 1D and 2D prescribed flow cases. That is, would the buffer super-droplets be subject to displacement like the rest? Or will the user be able to specify the position of buffer super-droplets once they "activate" them? I plan to look into this once I'm back from my summer break.

In terms of a test case from literature, the best reference I could think of is Prabhakaran et al. (2023). However, they use a LES and not a parcel model. @claresinger, I recall us discussing this paper and I wonder if you have a better reference in mind?

[slayoo]

Perhaps we could try with this paper as a source of an "example": https://doi.org/10.1175/1520-0450(1997)036%3C1449:CPTHSW%3E2.0.CO;2

[slayoo]

Probably this one could also be reproduced: https://doi.org/10.1098/rsta.2012.0086

[claresinger]

Both of these papers look like good examples to start with because the Parcel environment will be easier than a 1D case to start.

The LES is more complex than any of the environments we have in PySDM right now, but it will definitely an interesting paper to compare your results to in the future @jtbuch .

[claresinger]

Concrete next steps are outlined in issue #1387