open-atmos / PySDM

Pythonic particle-based (super-droplet) warm-rain/aqueous-chemistry cloud microphysics package with box, parcel & 1D/2D prescribed-flow examples in Python, Julia and Matlab
https://open-atmos.github.io/PySDM/
GNU General Public License v3.0
51 stars 28 forks source link

Time-dependent super-droplet sampling and modification #1331

Open jtbuch opened 2 months ago

jtbuch commented 2 months ago

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 commented 1 month ago

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 commented 1 month ago

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?