En este trabajo, consideramos un modelo de aguas someras de múltiples capas con densidad variable. Consiste en un sistema de ecuaciones hiperbólicas con productos no conservativos que tiene en cuenta las variaciones de presión debido a las fluctuaciones de densidad en un fluido estratificado. Se desarrolla un método de volumen finito de segundo orden que combina una técnica de reconstrucción hidrostática con un operador de reconstrucción de segundo orden MUSCL. El esquema está bien equilibrado para las soluciones estacionarias del lago en reposo. Además, se proporcionan indicaciones sobre cómo preservar una clase general de soluciones estacionarias correspondientes a un perfil de densidad estratificada. Se presentan algunos resultados numéricos, incluida la validación con datos de laboratorio que muestran la eficiencia y precisión del enfoque introducido aquí. Finalmente, se presenta una comparación entre dos estrategias de paralelización diferentes en GPU.