Reduced Order Models (ROMs) are simplified representations of large systems designed to reduce computational cost while retaining the essential features. Their key challenges are the identification of such features and the quantification of the error committed. Metric ROMs capture the geometric structure of the problem and provide geometry-aware models by identifying significant topological features of the nonlinear data manifold at disposal, addressing relevant limitations of traditional ROMs. This work proposes a nonlinear metric ROM within the Wasserstein metric space P2(Ω;W2) for gradient flow-based dynamics. Such a setting is particularly well-suited for handling both advective and diffusive behaviors, whose coexistence is notoriously problematic for the identification of accurate ROMs. Considering ρ0 ∈ P2(Ω) a probability density and vt ∈ L2(ρt;Ω), ∀t ∈ [0,T], a smooth vector field s.t. ρtvt ·ν ∂Ω = 0, the gradient flow is expressed with ρt ∈ P2(Ω) solving the continuity equation ∂tρt +∇·(ρtvt) = 0 in the distributional sense. Therefore, the parametric diffusion-advection-reaction PDE with Neumann boundary conditions is represented as a gradient flow with velocity vt(x) = −∇(log(ρ(x;p))+V(x;p)). We build a dynamical low-rank approximation of ρt(x;p) including its geometrical constraints, providing insights on the accuracy of the approximation. Since, commonly, nonlinear ROMs are built by applying linear reduction on a transformed manifold that is the result of a nonlinear transformation of the original snapshots, we propose a nonlinear ROM in which linear reduction is performed on an appropriate data manifold, found by exploiting the geometric assumptions of its gradient flow structure. We demonstrate that this improves the generalization capabilities in the entire parameter space, and apply the method to different scenarios.