Mixed-phase clouds, in which supercooled droplets and ice crystals coexist, remain one of the largest sources of uncertainty in climate prediction. One process that governs their glaciation is the Wegener–Bergeron–Findeisen (WBF) process, in which ice particles grow at the expense of evaporating droplets. This project studies WBF process with direct numerical simulation (DNS) using Lagrangian particle tracking. The solver is a pseudo-spectral code in which droplets and ice particles are tracked individually and coupled in both directions to the supersaturation field. Under the current NAISS CPU allocation (naiss2026-3-574-cpu), held with my main PhD supervisor Gaetano Sardina (Department of Mechanical Engineering, Chalmers University of Technology), we have established how turbulence, liquid–ice segregation and supersaturation fluctuations influence the glaciation time, assuming WBF as the main driver of this mechanism. The next stage of this project adds collision and riming — direct contact between two ice particles, and between ice particles and droplets — on top of the existing WBF physics. This model requires tracking many particles and adds a neighbor search over all particle pairs at every timestep, which makes the CPU version impractical for the planned cases. Due to this reason, the Arrhenius GH200 nodes are strongly needed for this next work. Even the largest planned configuration, 384³ grid resolution, needs only about 19 GB of field storage against 96 GB of high-bandwidth memory per GPU, so a full simulation fits on a single device with ample room left for large particle populations. The global FFT transposes that dominate communication cost in the distributed CPU version are therefore eliminated rather than merely accelerated. The particle work (e.g., interpolation, evaporation, deposition, and collision detection) is highly data-parallel and maps naturally onto the GPU. The GPU allocation will be used to port and validate the solver against the existing CPU code, implement and verify the collision and riming model, as well as run production cases under different initial conditions, Damköhler number, gravitational settling, initial droplet radius, and domain size, at resolutions up to 384³.