Numerical Modeling of Salinity Transport in the Panama Canal Neopanamax Locks Using DualSPHysics©

1. Introduction

The Panama Canal is one of the world's most significant hydraulic infrastructures, connecting the Atlantic and Pacific Oceans and facilitating a substantial portion of global maritime trade. In order to increase transit capacity and accommodate the growth of international shipping, the New Panama Canal locks, known as the Neopanamax locks, were constructed. These locks incorporate larger chambers and advanced water management systems designed to improve operational efficiency [1].

The operation of the Neopanamax locks has introduced new challenges related to water resources management and water quality. Among these, saltwater intrusion into Gatun Lake has become a growing concern because the lake serves as the primary freshwater source for canal operations and provides potable water to a substantial portion of the Panamanian population [2].

During lockage operations, water exchange between the lock chambers and the ocean results in the progressive transport of saline water into the lock system, generating density gradients that induce gravity currents and complex mixing processes between freshwater and seawater.

In the old locks of the Panama Canal, salinity levels gradually decrease from the ocean to Gatun Lake through the locks. Figure 1 shows the maximum and minimum salinity values of the Atlantic Ocean and the original Gatun Locks of the Panama Canal. The most critical phase of saltwater intrusion through the locks occurs when vessels are heading toward the ocean (down-lockages). As the vessel enters the locks, it displaces the satier water in the opposite direction of its movement.

Studies conducted during the design phase of the Neopanamax locks examined the potential for saltwater intrusion under various operating configurations. The results showed that salinity levels in Gatun Lake would remain below the limit of 0.5 ppt for brackish water [3,4,5,6,7,8]. Based on these results, the Panama Canal Master Plan concluded that implementing the Neopanamax locks would not significantly impact the lake’s water quality [1]. However, measurements obtained during the operational phase have revealed salinity increases exceeding the originally projected values, reaching over 1.0 ppt in some areas during the 2024 drought [9,10]. These observations have increased interest in understanding the physical mechanisms governing salinity transport within the lock system.

Several mitigation strategies were evaluated during the design process, including lock chamber flushing systems, deep pits for capturing saltwater and pneumatic barriers [7,8]. More recently, the Panama Canal Authority (ACP) has implemented operational strategies involving the use of the locks' water-saving basins to create temporary hydraulic barriers [10]. This reduces the advance of the saline front during specific lockage operations. Despite these measures, a more comprehensive understanding of the hydrodynamic processes controlling saltwater intrusion and its propagation towards Gatun Lake is still necessary.

Numerical modeling has become an essential tool for investigating density-driven mass transport phenomena. Among the available approaches, the Smoothed Particle Hydrodynamics (SPH) method has demonstrated a high capability for representing multiphase flows, highly deformable free surfaces, and complex interactions between fluids with different physical properties [11,12,13,14,15,16]. Due to its mesh-free Lagrangian nature, SPH is particularly suitable for simulating gravity currents, stratification processes, and mixing phenomena associated with saltwater intrusion.

In particular, the open-source solver DualSPHysics©, based on the weakly compressible SPH formulation of the Navier–Stokes equations, has become one of the most widely adopted tools in coastal and hydraulic engineering. It has been successfully applied to a broad range of coastal engineering problems, including wave propagation and breaking [17,18,19,20,21,22,23,24], hydraulic jumps [25,26], wave–structure and jet–wave interactions [27,28,29,30], as well as multiphase and density-driven flows [31,32]. Recent developments have significantly expanded the capability of SPH to represent complex hydraulic phenomena, including buoyancy-driven circulation, multi-component mixing, and density-driven flows. These advances have made SPH a suitable and robust tool for investigating salinity transport and stratification processes in large-scale hydraulic infrastructures [33].

The objective of this study is to investigate salinity transport processes within the Neopanamax locks of the Panama Canal through numerical simulations performed with the SPH-based solver DualSPHysics©. A two-dimensional numerical model representative of the Atlantic-side lock system was developed using projected salinity conditions [8] and considering different lock operation scenarios involving water exchange between adjacent chambers and between the locks and the Atlantic Ocean. Given the strategic importance of the Panama Canal for global maritime trade and the increasing pressure exerted by climate change on regional freshwater resources, this study provides insights into the physical mechanisms governing saltwater intrusion under future environmental conditions. The findings contribute to improving the understanding of density-driven transport processes in large-scale hydraulic infrastructures and provide a scientific basis for the development of effective operational and mitigation strategies aimed at preserving the freshwater resources of Gatun Lake while ensuring the long-term sustainability of canal operations.

2. Materials and Methods

2.1. Numerical Method

The numerical simulations were performed using the SPH method implemented in the open-source software DualSPHysics©. SPH is a mesh-free Lagrangian approach in which the fluid domain is discretized into particles that carry physical properties such as mass, density, velocity, and pressure. The governing equations of fluid motion are solved through kernel interpolation functions, allowing the simulation of highly deformable free surfaces and complex multiphase interactions. The general description of the method is omitted here, and the reader can refer to specialized textbooks [34]

The motion is represented by the Navier-Stokes Equations for a weakly compressible fluid. In the SPH Lagrangian formulation, the continuity and momentum equations for particle a, located at position r_a, are given by:

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_a}{\mathit{dt}}}={\rho }_a{\displaystyle \sum _b}{\displaystyle \frac{m_b}{{\rho }_b}}\left({\mathit{v}}_a-{\mathit{v}}_b\right)\bullet {\nabla }_a{W}_{\mathit{ab}}+D_a\\ {\displaystyle \frac{d{\mathit{v}}_a}{\mathit{dt}}}=-{\displaystyle \sum _b}{m}_b\left({\displaystyle \frac{P_a+P_b}{{\rho }_a{\rho }_b}}\right)\nabla W_{\mathit{ab}}+{\mathsf{\Gamma }}_a+g\end{array}$$

(1)

where ρ is the density, v is the velocity vector, P is the pressure, D_a represents a numerical diffusive term, Γ_adenotes the dissipation terms, and g is the gravity acceleration vector. Furthermore, W_ab = W(r_a − r_b , h) and (∙)_ab = (∙)_a − (∙)_b. The summations in Eq. (1) are extended to all the particles b located within the circular domain centered on a and of radius 2h (with h the smoothing length), where the kernel function W_ab is defined. Here, the adopted kernel function W_ab is the cubic-spline kernel function [35].

The simulations were performed using the multiphase non-Newtonian formulation implemented in DualSPHysics©, based on the Herschel–Bulkley–Papanastasiou (HBP) constitutive model [36]. In a non-Newtonian fluid, the shear stresses, τ, are given by:

$$\mathit{\tau }{=2\eta }_{\mathit{app}}\dot{\mathit{\gamma }}$$

(2)

where η_app is the apparent viscosity and $\dot{\mathit{\gamma }}$is the symmetric strain rate tensor:

$$\dot{\mathit{\gamma }}$$

$$={\displaystyle \frac{1}{2}}\left[\nabla \mathit{v}+{\nabla }^{\mathit{T}}\mathit{v}\right]-{\displaystyle \frac{1}{3}}\left(\nabla .\mathit{v}\right)\mathit{I}$$

(3)

According to the HBP constitutive equation used in DualSPHysics©, η_app is expressed as:

$${\eta }_{\mathit{app}}=\eta {\left|\dot{\gamma }\right|}^{n_{\mathit{HB}}-1}+{\displaystyle \frac{{\tau }_y}{\left|\dot{\gamma }\right|}}\left[1-e^{-m_p\left|\dot{\gamma }\right|}\right]$$

(4)

In this formula, $\left|\dot{\gamma }\right|$ represents the magnitude of the symmetric strain rate tensor, τ_y denotes the yield strength of the material, and n_HB and m_p represent the power law index and exponential shear rate growth parameter, respectively. This generalized rheological formulation is capable of representing different time-independent non-Newtonian behaviors through appropriate values of n_HB, m_p, and τ_y.

In the present study, freshwater and saline water were assumed to behave as Newtonian fluids. Accordingly, the rheological parameters were set to n_HB = 1, m_p = 0, and τ_y = 0, so that the HBP model reduces to the Newtonian constitutive law, $\mathit{\tau }=2\eta \dot{\mathit{\gamma }}.$ The adopted multiphase formulation therefore allows two interacting Newtonian fluids with different densities to be represented, enabling the simulation of density-driven transport, stratification, and mixing processes associated with saltwater intrusion.

2.2. Salinity Transport Model Validation

The numerical salinity transport model was validated against the classical lock-exchange experiment reported by Shin et al. [37], which has been widely used as a benchmark for density-driven gravity currents. The experiment consists of the sudden release of two fluids with different densities initially separated by a vertical gate in a rectangular channel. After gate removal, the denser saline water propagates along the bottom while the lighter freshwater moves above it, producing a bidirectional gravity current.

The validation case was reproduced using the Multi-phase Non-Newtonian solver available in DualSPHysics©. The computational domain and fluid properties were defined according to the experimental configuration reported by Shin et al. [37]. The experimental configuration consisted of a rectangular channel with length L = 2 m and depth H = 0.2 m. Freshwater and saline water were assigned densities of ρ₁ = 1000 kg/m³ and ρ₂ = 1007 k g/m³, respectively, corresponding to a density ratio of γ = ρ₁/ρ₂ = 0.993. This configuration reproduces the stratification and mixing processes characteristic of lock-exchange flows. To facilitate comparison, the results were expressed in terms of the dimensionless front position x^* = x/H, where x is the front position. Following the original study, dimensionless time was expressed as $t^{*}=t\sqrt{g(1-\mathsf{\gamma }/\mathrm{H}})$, where t is time and g is the acceleration due to gravity.

The salinity transport model used for the validation phase was constructed in DalSPHysics©, with an inter-particle distance of 0.001 m. This configuration resulted in 785,000 particles, with a computing time of 10 hours on a computer equipped with an NVIDIA© GeForce RTX 2060 GPU.

Figure 2 compares the temporal evolution of the density interface obtained from the SPH simulation with the experimental observations. The numerical model successfully reproduces the development of the gravity current, including the propagation of the saline front and the evolution of the density interface.

Figure 3 presents the comparison between the dimensionless gravity current front position obtained experimentally and that predicted by the SPH model. The numerical results closely follow the experimental measurements throughout the entire propagation process, accurately reproducing both the propagation rate and the overall evolution of the gravity current front. Only minor deviations are observed at a few dimensionless times, while the overall agreement remains excellent. This qualitative agreement is further supported by the statistical indicators, yielding a Root Mean Square Error (RMSE) of 0.0738, a coefficient of determination (R²) of 0.9994, and a Mean Absolute Percentage Error (MAPE) of 8.64%. These results confirm that the SPH model accurately reproduces the lock-exchange experiment and provides a reliable representation of density-driven gravity currents, supporting its application to salinity transport simulations in the Panama Canal locks.

2.3. Numerical Models of the Simulated Scenarios

To simulate salinity transport in the Panama Canal locks using DualSPHysics©, the dimensions of the Post-Panamax locks on the Atlantic side were used as a reference (Figure 4), as presented by the Panama Canal Authority [8]. Using these dimensions ensured that the numerical simulations accurately reflected the hydraulic conditions and the dynamics of salinity transport in the Panama Canal locks, allowing for a more precise evaluation of the behavior of gravity currents and the mixing processes between saltwater and freshwater. The design included a complete representation of the lock chamber structures, the operational water levels, and the hydraulic connections to the Atlantic Ocean and Gatun Lake. The numerical simulations were performed in two dimensions (2D), which significantly reduced the computational load without compromising the accuracy of the physical phenomena represented.

A scale of 1:500 was applied to the horizontal axis and 1:50 to the vertical axis, with the aim of improving the visualization of the salinity transport below the vessels, as well as the vertical mixing behind them. Additionally, a waiting time of 5 min was implemented while the vessels transited through the lock chambers, in order to correspond with the gate-opening process. The vessels then moved at a speed of 2 knots, slowing down until they came to a halt at the end of the chamber. This illustrates the movement of the vessels and aims to realistically simulate the hydrodynamic effect generated by their forward motion on the mixing pattern and propagation of the salinity front within the chambers.

The DualSPHysics© salinity transport models of the simulated scenarios used an inter-particle distance of 0.002 m, resulting in 1,200,000 to 2,000,000 particles, with a computing time of 2 to 3 days on a computer equipped with an NVIDIA© GeForce RTX 2060 GPU. For the simulations, the maximum average salinity concentrations shown in Figure 4 were used as a baseline. These concentrations correspond to the low water levels in the Neopanamax lock chambers on the Atlantic side and are based on projections by the Panama Canal Authority, as presented in a WL Delft Hydraulics report [8]. These values correspond to the salt concentration range in parts per thousand (ppt) for the lock chambers at the end of the dry season in 2060. The specific values used in this context are from Scenario 1, in which the water saving basins of the new Neopanamax locks operate at full capacity and water losses caused by lock operations are offset by additional water supplies to Gatun Lake, in addition to current supplies [8]. The high water levels in the locks were reconstructed from the baseline low water levels of Figure 4 through mass conservation analysis of the locks equalization processes.

Five numerical scenarios were designed to investigate the transport and mixing of saline water under the Panama Canal Neopanamax locks' representative operating conditions. All simulations were performed using the validated Multi-phase Non-Newtonian solver implemented in DualSPHysics©. These scenarios reproduce different lock operations involving adjacent chambers, the Atlantic Ocean and Gatun Lake, enabling the impact of density-driven currents and vessel motion on salinity intrusion to be assessed.

The first two scenarios (Case 1 and Case 2) simulate the opening of gates between adjacent lock chambers with different salinity levels, and the vessels transit through the locks. These cases assess the progression of the density interface and the mixing process triggered by salinity gradient and the vessels' transit between adjacent chambers. The third and fourth scenarios represent the interaction between the Atlantic Ocean and the outermost lock chamber during vessel departure (Case 3) and arrival (Case 4), respectively. These simulations investigate the influence of vessel direction on the propagation of saline intrusion from the ocean into the lock system under realistic operating conditions. Scenarios fifth and sixth evaluate the transport of salinity between the lock chamber adjacent to Gatun Lake and the freshwater reservoir (Case 5), as well as a mitigation strategy consisting of a trapezoidal retention pit (Case 6). The aim in Case 6 is to determine the effectiveness of the pit in retaining the denser saline water and reducing saltwater intrusion towards the lake.

All scenarios were analyzed in terms of the temporal evolution of the salinity field, density interface and gravity-current propagation. Particular attention was given to identifying the transport mechanisms and evaluating the effectiveness of the proposed mitigation strategy.

3. Results and Discussion

3.1. Salinity Transport Between Adjacent Lock Chambers

Case 1 and Case 2 represents the opening of the gates between adjacent lock chambers and the vessels transit through the locks M and L (Case 1), and L and K (Case 2).

The chambers initially contain water with different salinity levels, generating a density gradient that drives the gravity current after gate opening. In Case 1, the salinity S and density ρ values in both chambers at the start of the simulation were S = 12 ppt (ρ = 1007.44 kg/m³) for lock M and S = 21 ppt (ρ = 1014.29 kg/m³) for lock L. In Case 2 the salinity and density values at the start of the simulation were S = 21 ppt (ρ = 1014.29 kg/m³) for lock L and S = 27 ppt (ρ = 1018.84 kg/m³) for lock K. The schematic configuration for Case 1 is shown in Figure 5.

Figure 6a shows the initial state before the gates are opened. Lock M contains lower-salinity water, shown in blue, while lock L contains higher-salinity water, shown in red. Immediately after the gates begin to open, the denser saline water moves along the bottom of the chamber while the lighter freshwater moves above it, creating a classic lock-exchange gravity current (Figure 6b). After the 5 min gate-opening period has elapsed, the vessel (blue rectangle) begins to move toward the adjacent chamber, while the saltier water has flowed beneath the vessel along almost its entire length (Figure 6c). As the vessel moves through the adjacent lock, the flow of saltwater beneath the vessel reaches the opposite wall of the upper lock and begins to move upward (Figure 6d). Finally, the vessel enters the lock and slows down until it comes to a complete stop in the lower lock chamber, while, behind it, the salinity concentration is distributed fairly evenly throughout the entire volume of water (Figures 6e–6g). In Case 2, the initial salinity difference is slightly different from that of Case 1, but the physical transport mechanism is very similar to that of Case 1.

3.2. Salinity Intrusion During Lock-Ocean and Ocean–Lock Vessel Transits

Case 3 simulates the opening of the gate between the Atlantic Ocean, O, and lock K, followed by the vessel departing the lock towards the ocean. In this case, the salinity and density values at the start of the simulation were S = 27 ppt (ρ = 1018.84 kg/m³) for lock K and S = 32 ppt (ρ = 1022.68 kg/m³) for the Atlantic Ocean, O. This operating condition represents the main pathway for saltwater to intrude into the lock system. Figure 7 illustrates the schematic configuration for this scenario.

Figure 8a shows the initial stage of the simulation. Lock K containing lower-salinity water (shown in blue) and the Atlantic Ocean containing higher-salinity water (shown in red), are clearly visible. Immediately after the gates begin to open, the denser ocean water enters the lock chamber as a bottom gravity current and moves progressively towards the interior of the chamber (Figures 8b–8c). Over time, turbulent mixing, aided by vessel movement, redistributes the saline water until an almost homogeneous salinity field is achieved (Figures 8d–8h). The results suggest that vessel departure favors relatively large saltwater intrusion into the lock chamber.

Case 4 analyses the opposite operating condition to Case 3, where the vessel enters Lock K from the Atlantic Ocean. At the beginning of the gate opening, the vessel occupies the transition zone between the ocean and Lock K (Figure 9a). The saline front then propagates along the chamber bottom, creating a return current of lighter water that encounters the bow of the vessel (Figures 9b–9d). After the 5 min gate opening period, the vessel's movement pushes the lighter water below and behind it (Figures 9e–9f). The vessel's presence partially obstructs the flow, delaying the advance of the dense current and reducing the volume of saline water entering the chamber (Figure 9g). Although complete mixing eventually occurs, the intrusion process is noticeably slower than in Case 3, demonstrating the influence of the vessel's position on density-driven transport.

3.3. Salinity Intrusion During Lake-Lock Transits and Performance of Retention Pit

The main mechanism by which salinity enters Gatun Lake, the down-lockages (Case 5), is illustrated in Figure 10. In this case, the salinity and density values at the start of the simulation were S = 12 ppt (ρ = 1000 kg/m³) for Gatun Lake, G, and S = 32 ppt (ρ = 1007.44 kg/m³) for Lock M.

At the start of the gate opening process, the vessel occupies the transition zone between Gatun Lake and Lock M (see Figure 10a). The saline front then propagates along the lake bottom (Figure 10b). After the 5 min gate opening period, the vessel's movement pushes the saline water below and behind it (Figures 10c-10d). Eventually, most of the saltier water in Lock M is transferred to Gatun Lake (Figure 10e).

Case 6 investigates the effectiveness of a trapezoidal retention pit installed between Lock M and Gatun Lake, the freshwater reservoir. The pit is intended to retain the denser saline water before it reaches the lake, thereby reducing saltwater intrusion. The conceptual configuration of this mitigation strategy is shown in Figure 11.

Figure 12 shows the temporal evolution of salinity transport using a retention pit. At the beginning of the gate opening, the vessel occupies the zone above the trapezoidal retention pit (Figure 12a). Once the gate opening process is initiated, and beyond, the saline water descends into the pit under the action of gravity (Figures 12b–12e). It eventually reaches and accumulates at the far end of the pit (Figures 12f–12g). It then returns to the beginning of the pit (Figure 12h), and finally becomes concentrated in the bottom depression (Figure 12i). Only a small amount of the saline water continues towards Gatun Lake, while most of the denser water is retained inside the pit. This behavior substantially reduces the propagation of saltwater into the freshwater reservoir, demonstrating the possible effectiveness of the proposed passive hydraulic structure as a salinity mitigation measure.

It is important to note that, due to the use of different scales in the horizontal and vertical directions, the turbulence mechanism is necessarily distorted; therefore, especially in this case of the retention pit, the results should be regarded primarily as qualitative. To obtain a more accurate result, it is essential to use the same numerical scale in both directions. This may require a larger number of particles and the use of computers with multiple GPUs. Alongside 3D simulations, this would undoubtedly be the next stage in saltwater propagation studies.

4. Conclusions

This study demonstrates the capability of the SPH method implemented in DualSPHysics® to simulate density-driven salinity transport in the Panama Canal lock system. The numerical model successfully reproduced gravity-current dynamics and salinity mixing processes under representative lock operating conditions, providing a robust tool for analyzing complex free-surface flows involving fluids with different densities.

The simulations of adjacent lock chambers showed that gate opening generates rapid gravity-current propagation and efficient nearly homogeneous mixing between water masses of different salinity.

The simulations involving the Atlantic Ocean demonstrated that vessel direction significantly influences saltwater intrusion into the lock system. When the vessel exits toward the ocean, a larger volume of saline water enters the lock chamber, whereas a vessel entering from the ocean partially obstructs the density current, reducing and delaying saltwater intrusion. These findings indicate that vessel position should be considered when evaluating salinity transport during lock operations.

The proposed trapezoidal retention pit proved effective in mitigating saltwater intrusion toward Gatun Lake. The denser saline water remained preferentially trapped within the pit, substantially reducing the downstream propagation of the saline front. This passive hydraulic solution represents a promising alternative for limiting salinity intrusion without modifying normal lock operations.

Overall, the results demonstrate that SPH numerical simulations provide an effective framework for evaluating hydraulic mitigation strategies and improving the understanding of salinity transport processes in navigation lock systems. The methodology developed in this study can support future investigations aimed at optimizing lock operations and designing infrastructures to reduce saltwater intrusion into freshwater reservoirs. The next phase of the study should involve using the same numerical scale in both horizontal and vertical directions, as well as 3D simulations utilizing computers with multiple GPUs.