Solver recommendation for large advection-diffusion (500k cells)

Hi,
I'm at the New University of Lisbon (NOVA), working on a collaborative project with Dalhousie University and the Canadian DFO.
Background: I've successfully used CVODE for a diagenetic model (stiff biogeochemistry + vertical advection-diffusion in sediment layers). It works brilliantly - thank you to the development team!
New problem: I'm extending a related application that simulates dissolved substance transport (e.g. ammonia) from fish farms. The spatial domain is large:

100-120 × 100-120 XY cells
20-30 vertical (Z) layers
Up to ~500,000 cells total

The physics is advection-diffusion between cells - not stiff, just large. Currently using explicit Euler/RK but hitting stability limits with diffusion.
Question: Which SUNDIALS solver would you recommend for this scale?

Memory efficiency is a key concern given the matrix size.
Thanks for any guidance,
João G. Ferreira

Start with the SUNDIALS solver overview and the CVODE documentation, focusing on the stated 500,000-cell advection-diffusion case and its memory constraints. Compare the relevant solver guidance against the reported explicit Euler/RK stability problem. Done means documenting a justified solver recommendation for this scale.

Written by the indexing model from the issue text.