diff --git a/PySDM/formulae.py b/PySDM/formulae.py
index 853000eb7..c86763832 100644
--- a/PySDM/formulae.py
+++ b/PySDM/formulae.py
@@ -54,6 +54,8 @@ def __init__(  # pylint: disable=too-many-locals
         optical_albedo: str = "Null",
         optical_depth: str = "Null",
         particle_shape_and_density: str = "LiquidSpheres",
+        liquid_particle_shape: str = "LiquidSpheres",
+        ice_particle_shape: str = "IceSpheres",
         terminal_velocity: str = "GunnKinzer1949",
         air_dynamic_viscosity: str = "ZografosEtAl1987",
         bulk_phase_partitioning: str = "Null",
diff --git a/PySDM/physics/particle_shape_and_density/liquid_spheres.py b/PySDM/physics/particle_shape_and_density/liquid_spheres.py
index 8b31b5682..61eee9d57 100644
--- a/PySDM/physics/particle_shape_and_density/liquid_spheres.py
+++ b/PySDM/physics/particle_shape_and_density/liquid_spheres.py
@@ -26,6 +26,3 @@ def volume_to_mass(const, volume):
     def radius_to_mass(const, radius):
         return const.rho_w * const.PI_4_3 * np.power(radius, const.THREE)
 
-    @staticmethod
-    def reynolds_number(_, radius, velocity_wrt_air, dynamic_viscosity, density):
-        return 2 * radius * velocity_wrt_air * density / dynamic_viscosity
diff --git a/PySDM/physics/particle_shape_ice/__init__.py b/PySDM/physics/particle_shape_ice/__init__.py
new file mode 100644
index 000000000..53f1e2896
--- /dev/null
+++ b/PySDM/physics/particle_shape_ice/__init__.py
@@ -0,0 +1,5 @@
+"""
+particle shape and density relationships for ice particles
+"""
+
+from .ice_spheres import IceSpheres
diff --git a/PySDM/physics/particle_shape_ice/columnar.py b/PySDM/physics/particle_shape_ice/columnar.py
new file mode 100644
index 000000000..31a5b37c9
--- /dev/null
+++ b/PySDM/physics/particle_shape_ice/columnar.py
@@ -0,0 +1,12 @@
+"""
+columnar particles with constant density of ice
+"""
+
+class Columnar:
+    def __init__(self, const, _):
+        self.mass_density = const.rho_i
+
+    @staticmethod
+    def supports_mixed_phase(_=None):
+        return True
+
diff --git a/PySDM/physics/particle_shape_ice/ice_spheres.py b/PySDM/physics/particle_shape_ice/ice_spheres.py
new file mode 100644
index 000000000..b7ecc450c
--- /dev/null
+++ b/PySDM/physics/particle_shape_ice/ice_spheres.py
@@ -0,0 +1,63 @@
+"""
+spherical particles with constant density of ice
+"""
+from PySDM.physics.trivia import Trivia
+
+class IceSpheres:
+    def __init__(self, const, _):
+        pass
+
+    @staticmethod
+    def supports_mixed_phase(_=None):
+        return True
+
+    @staticmethod
+    def mass_to_volume(const, mass):
+        return (
+            Trivia.volume_of_density_mass(const, const.rho_i, mass)
+        )
+
+    @staticmethod
+    def volume_to_mass(const, volume):
+        return (
+            Trivia.mass_of_density_volume(const.rho_i, volume)
+        )
+
+    @staticmethod
+    def volume_to_radius(const, volume):
+        return (
+            Trivia.radius(const, volume)
+        )
+
+    @staticmethod
+    def radius_to_volume(const, radius):
+        return (
+            Trivia.sphere_radius_to_volume(const, radius)
+        )
+
+    @staticmethod
+    def radius_to_mass(const, radius, mass_density):
+        return (
+            Trivia.sphere_radius_to_volume(const, radius) * const.rho_i
+        )
+
+    @staticmethod
+    def radius_to_area(const, radius):
+        return(
+            Trivia.area(const, radius)
+        )
+
+    @staticmethod
+    def mass_to_capacity(const, mass):
+        return(
+            Trivia.sphere_mass_to_radius(const, mass, const.rho_i,)
+        )
+
+    @staticmethod
+    def maximum_diameter(const, mass):
+        return(
+            Trivia.sphere_mass_to_radius(const, mass, const.rho_i) * 2
+        )
+
+
+
diff --git a/PySDM/physics/particle_shape_liquid/__init__.py b/PySDM/physics/particle_shape_liquid/__init__.py
new file mode 100644
index 000000000..a35c5f0c4
--- /dev/null
+++ b/PySDM/physics/particle_shape_liquid/__init__.py
@@ -0,0 +1,6 @@
+"""
+particle shape and density relationships for liquid particles
+"""
+
+from .liquid_spheres import LiquidSpheres
+from .porous_spheroids import PorousSpheroid
diff --git a/PySDM/physics/particle_shape_liquid/liquid_spheres.py b/PySDM/physics/particle_shape_liquid/liquid_spheres.py
new file mode 100644
index 000000000..3351bc253
--- /dev/null
+++ b/PySDM/physics/particle_shape_liquid/liquid_spheres.py
@@ -0,0 +1,56 @@
+"""
+spherical particles with constant density of ice
+"""
+from PySDM.physics.trivia import Trivia
+
+class LiquidSpheres:
+    def __init__(self, const, _):
+        pass
+
+    @staticmethod
+    def mass_to_volume(const, mass):
+        return (
+            Trivia.volume_of_density_mass(const, const.rho_w, mass)
+        )
+
+    @staticmethod
+    def volume_to_mass(const, volume):
+        return (
+            Trivia.mass_of_density_volume(const.rho_w, volume)
+        )
+
+    @staticmethod
+    def volume_to_radius(const, volume):
+        return (
+            Trivia.radius(const, volume)
+        )
+
+    @staticmethod
+    def radius_to_volume(const, radius):
+        return (
+            Trivia.sphere_radius_to_volume(const, radius)
+        )
+
+    @staticmethod
+    def radius_to_mass(const, radius, mass_density):
+        return (
+                Trivia.sphere_radius_to_volume(const, radius) * const.rho_w
+        )
+
+    @staticmethod
+    def radius_to_area(const, radius):
+        return (
+            Trivia.area(const, radius)
+        )
+
+    @staticmethod
+    def mass_to_capacity(const, mass):
+        return (
+            Trivia.sphere_mass_to_radius(const, mass, const.rho_w, )
+        )
+
+    @staticmethod
+    def maximum_diameter(const, mass):
+        return (
+                Trivia.sphere_mass_to_radius(const, mass, const.rho_w) * 2
+        )
\ No newline at end of file
diff --git a/PySDM/physics/particle_shape_liquid/porous_spheroids.py b/PySDM/physics/particle_shape_liquid/porous_spheroids.py
new file mode 100644
index 000000000..94d6cf3d6
--- /dev/null
+++ b/PySDM/physics/particle_shape_liquid/porous_spheroids.py
@@ -0,0 +1,10 @@
+"""
+for mixed-phase microphysics as in
+[Shima et al. 2020](https://doi.org/10.5194/gmd-13-4107-2020)
+"""
+
+
+class PorousSpheroid:  # pylint: disable=too-few-public-methods
+    @staticmethod
+    def supports_mixed_phase(_=None):
+        return True
diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py
index 8bfb0f771..97852e3db 100644
--- a/PySDM/physics/trivia.py
+++ b/PySDM/physics/trivia.py
@@ -12,21 +12,36 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def volume_of_density_mass(rho, m):
-        return m / rho
+    def volume_of_density_mass(rho, mass):
+        return mass / rho
 
+    @staticmethod
+    def mass_of_density_volume(rho, volume):
+        return volume * rho
+
+    # TODO: change name to sphere_volume_to_radius
     @staticmethod
     def radius(const, volume):
         return np.power(volume / const.PI_4_3, const.ONE_THIRD)
 
+    # TODO: change name to sphere_radius_to_area
     @staticmethod
     def area(const, radius):
         return const.PI * const.FOUR * np.power(radius, const.TWO)
 
     @staticmethod
-    def volume(const, radius):
+    def sphere_radius_to_volume(const, radius):
         return const.PI_4_3 * np.power(radius, const.THREE)
 
+    @staticmethod
+    def sphere_radius_to_mass(const, radius, density):
+        return const.PI_4_3 * np.power(radius, const.THREE) * density
+
+    @staticmethod
+    def sphere_mass_to_radius(const, mass, density):
+        return const.PI_4_3 * np.power(radius, const.THREE) * density
+
+    # TODO: remove this. There is already sphere_radius_to_area
     @staticmethod
     def sphere_surface(const, diameter):
         return const.PI * diameter**2