wet radius initialisation (with examples from 1D kinematic example)

[jbarr444]

Hi Sylwester,

Working with the Shipway/Hill example. I find strange behavior with the particle size distribution at time 0, where the wet and dry spectra instantaneously differ dramatically at a certain height corresponding to the maximum initial relative humidity. I think this growth is far beyond equilibrium, so I'm puzzled. Can you explain what's going on under the hood?

wet spectrum: https://github.com/open-atmos/PySDM/assets/166780711/b306d37f-4e0b-471e-abe4-6d2616b5336b
RH vs z: https://github.com/open-atmos/PySDM/assets/166780711/7688e95e-56f1-4fbd-9f0f-39c72bb9afa7

Thanks,
Jason

[slayoo]

Hi Jason,

Thanks for reporting it. Could you elaborate on how the "dry spectra instantaneously differ"?

Trying with a much lower kappa could help to check if this is not a matter of high hygroscopicity?

The way the RH profile is plotted might be a bit misleading here, as we do not spatially interpolate RH for super particles; plotting a stepwise constant profile could help to connect the two plots.

HTH,
Sylwester

[jbarr444]

Hi Sylwester, thanks for the reply,

Here's the dry spectrum: https://github.com/open-atmos/PySDM/assets/166780711/bb6f3520-62f0-4291-8f24-afab0d212049

Notice how there's no 'kink' at 740 m. I'm wondering how the initial wet spectrum is calculated from the initial dry spectrum, as clearly it is not assumed that they are the same to begin with. I would think perhaps then you start at equilibrium with respect to condensation/evaporation, but this also does not seem to be the case, though I could be wrong about that.

Lowering the hygroscopicity from 0.9 to 0.1 shifts the wet sizes closer to the dry spectrum, but the 'kink' remains:
https://github.com/open-atmos/PySDM/assets/166780711/1a51676e-368a-4135-a731-9a8019d67367

[slayoo]

The dry spectrum is set to constant with altitude in this example, so the only variability will originate from random sampling - consistent with your plot.

The "kink" in the wet spectrum comes from exactly what you indicate - the equilibrium wrt ambient RH, which is solved for at initialization here:

• PySDM/examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py

  Lines 95 to 103 in 4de8327

  	self.attributes = self.env.init_attributes(
  	spatial_discretisation=spatial_sampling.Pseudorandom(),
  	spectral_discretisation=spectral_sampling.ConstantMultiplicity(
  	spectrum=settings.wet_radius_spectrum_per_mass_of_dry_air
  	),
  	kappa=settings.kappa,
  	collisions_only=not settings.enable_condensation,
  	z_part=settings.z_part,
  	)

• PySDM/PySDM/environments/kinematic_1d.py

  Lines 73 to 78 in 4de8327

  	r_wet = equilibrate_wet_radii(
  	r_dry=r_dry,
  	environment=self,
  	cell_id=attributes["cell id"],
  	kappa_times_dry_volume=attributes["kappa times dry volume"],
  	)

• PySDM/PySDM/initialisation/equilibrate_wet_radii.py

  Lines 16 to 25 in 4de8327

  	def equilibrate_wet_radii(
  	*,
  	r_dry: np.ndarray,
  	environment,
  	kappa_times_dry_volume: np.ndarray,
  	f_org: np.ndarray = None,
  	cell_id: np.ndarray = None,
  	rtol=default_rtol,
  	max_iters=default_max_iters,
  	):

HTH,
S.

[jbarr444]

@slayoo Hi, hoping to hear back from you about this. Thanks so much.

[slayoo]

Hi @jbarr444,

The state variable triplet used in PySDM is: (dry air potential temperature, vapour mixing ratio, dry air density).

RH is computed in

PySDM/PySDM/backends/impl_numba/methods/physics_methods.py

Lines 56 to 62 in bafd7f8

	T[i] = ff.state_variable_triplet__T(rhod[i], thd[i])
	p[i] = ff.state_variable_triplet__p(
	rhod[i], T[i], water_vapour_mixing_ratio[i]
	)
	RH[i] = ff.state_variable_triplet__pv(
	p[i], water_vapour_mixing_ratio[i]
	) / ff.saturation_vapour_pressure__pvs_Celsius(T[i] - ff.constants.T0)

with:

• T = T(rhod, thd)
• p = p(rhod, T, mixing_ratio)
• pv = pv(p, mixing_ratio)
• RH = pv/pvs(T)

from https://github.com/open-atmos/PySDM/blob/main/PySDM/physics/state_variable_triplet/libcloudphplusplus.py and https://github.com/open-atmos/PySDM/tree/main/PySDM/physics/saturation_vapour_pressure

HTH,
S.

[jbarr444]

Hi @slayoo

Thanks for the tip. I can see now that PySDM is doing it right, and I was wrong by not factoring in the effect of humidity on pressure/temperature. Surprising how much of an effect it had. It appears the Shipway and Hill paper might've been flawed for having an initial ambient supersaturation..?

Case closed as far as I'm concerned.

Thanks again,
Jason

[slayoo]

@jbarr444, the 1D kinematic setup of Shipway & Hill assumes (quoting the paper) that the "temperature field is kept fixed throughout the simulation", but it is not clarified which temperature it is. The KiD code (https://github.com/BShipway/KiD) also uses two distinct density profiles - different one for microphysics, and a different one for advection. In PySDM, we assumed that the constant-in-time temperature is the dry-air-potential temperature, and we use the same hydrostatic dry-air-density profile both for advection and for thermodynamics. Anyhow, with the microphysics decoupled from thermodynamics (constant-in-time temperature profile), the thermodynamics are by definition not physical, so no matter how one interprets it, it will hardly be "doing it right" :)

[jbarr444]

@slayoo Interesting. Thanks for the clarification!