  We derive a linearized rotating shallow water system modeling tides, which
can be discretized by mixed finite elements. Unlike previous models, this model
allows for multiple layers stratified by density. Like the single-layer
case~\cite{kirby2021preconditioning} a weighted-norm preconditioner gives a
(nearly) parameter-robust method for solving the resulting linear system at
each time step, but the all-to-all coupling between the layers in the model
poses a significant challenge to efficiency. Neglecting the inter-layer
coupling gives a preconditioner that degrades rapidly as the number of layers
increases. By a careful analysis of the matrix that couples the layers, we
derive a robust method that requires solving a reformulated system that only
involves coupling between adjacent layers. Numerical results obtained using
Firedrake confirm the theory.
