Abstract

This paper deals with the initial modeling of water salinity and its diffusion into the lakes during lock operation on the Panama Canal. A hybrid operational model was implemented using the AnyLogic software simulation environment. This was accomplished by generating an operational discrete-event simulation model and a continuous simulation model based on differential equations, which modeled the salinity diffusion in the lakes. This paper presents that unique application and includes the effective integration of lock operations and its impact on the environment.

1. Introduction

The Panama Canal currently carries 4 percent of the world’s trade goods, and it is an important competitor in some very important shipping routes. For example, the canal currently handles about 16% of the United States maritime trade, and more than 25% of the containerized trade between north east Asia and the east coast of the United States [1, 2]. Within the Republic of Panama, the Canal is responsible for the growth of the terminal cities of Panama and Colon (these cities and their metropolitan areas have 70% of the population of the Republic of Panama).

The navigational channel of the Panama Canal is about 50 miles (80 km) long and extends from the Caribbean Sea to the Pacific ocean (see Figure 1). The canal includes Gatun Lake, a fresh water lake which has a surface area of 436 Km2 (168.4 mi2) at 26 m (85 ft) above mean sea level. A series of locks on the waterway is used to raise and lower ships transiting to and from Gatun Lake to the oceans by moving water into and out of the lake using the force of gravity [3, 4]. In addition, Gatun Lake is utilized to provide fresh water to the terminal cities of Panama and Colon.

Currently the canal is operating at almost 95% of its maximum sustainable capacity which limits its ability to capture the increasing demand that will occur within the next 50 years [5]. Additionally, port operators at both sides of the Canal are betting on larger and more efficient post-Panamax ships (e.g., “ships that do not fit in the canal, such as supertankers and the largest modern container ships”) [6]. Therefore, Panama has started the expansion of the Panama canal. This project will double the capacity of the Panama Canal by 2014 by allowing more and larger ships to transit.

The Panama Canal expansion through a third set of locks presents new challenges. One challenge is the potential increase in the salinity of Gatun Lake above permissible levels. Using this scenario we designed and propose to implement a system that allows us to estimate increases and/or decreases in salinity in Gatun Lake so that we can alternatively design and implement plans to reduce or mitigate adverse effects.

To understand this operational system and its complexities, we have to both describe and model its structure. Therefore, we have created a basic hybrid simulation model of the current Panama Canal operations that also shows the diffusion of the saline ocean water into Gatun Lake. The behavior of the Panama Canal, taking into consideration the transit of ships and the diffusion of salinity, is determined by the interaction of discrete event and continuous dynamics. We decided to combine explicitly a discrete event model (which models the operations and the transit of ships) and a continuous model which includes the respective differential equations to model the exchange of water and salt based on the mass balance and salt balance dynamics. AnyLogic provides the mechanisms to model this type of hybrid systems [7]. AnyLogic allows us, using hybrids statecharts (a form of Hybrid Automata), to combine, explicitly, the basic model of differential equations with the basic model of discrete event systems, which are finite state automata. The objective is to have a simulation model which can model the current transit rates of the Panama Canal (based on appropriate 95% confidence intervals) and the salinity levels with a mean absolute error lower than 0.02 ppt as recommended by Subject Matter Experts (SMEs).

This paper discusses the preliminary phase of this project where the current conditions of the Panama Canal are modeled using a basic hybrid model in AnyLogic. A potential increase in salinity due to the expansion of the Panama Canal with a new set of locks will be considered in a future paper.

2. Salinity of the Lakes

Salinity refers to the mass quantity of dissolved salts per unit of water mass or water volume (1 unit = 1 liter). Seawater’s salinity (𝑆) amounts to 35 parts per thousand (ppt). The chloridity (Cl), which is sometimes used, represents the mass quantity of chloride ions per unit of water mass or water volume. The Cl of fresh water should not exceed 0.2 to 0.25 ppt. This fresh water limit corresponds to a salinity value of 0.4 to 0.5 ppt. The relationship between 𝑆 and Cl is: 𝑆=0.03+1.905*Cl (valid for 𝑆 and Cl > 1 ppt) [8].

The salinity diffusion was modeled using the exchange of mass transfer [9]. This involved the study of different volumes and salinity gradients of the Panama Canal System: Water systems and Locks. In addition, it was complemented by a data collection using historical data provided by the Panama Canal Authority (and data collected by the team). For convenience reference, the six locks (see Figure 1) are identified by the following notation [10]: L1: lowest lock on the Pacific side;L2: lock between Miraflores Lake and L1;L3: lock between Miraflores Lake and Gatun Lake (Pedro Miguel);L4: highest lock directly connected to Gatun Lake;L5: middle lock on Atlantic side;L6: lowest lock connected to Atlantic ocean.

The Panama Canal has two lakes: Gatun and Miraflores. Water from these two lakes is used for the Panama Canal System to fill the navigation locks. Salt water from the Pacific and Atlantic Oceans gets added to the lakes during the transit of the ships. In addition, water from the lakes is lost to the sea during the same process. The Gatun Lake supplies fresh water to the population of Panama and Colon Cities for drinking purposes. The Miraflores Lake has a level of salinity which is already considered “brackish” water (i.e., brackish water is water that has more salinity than fresh water, but not as much as seawater) [8, 9, 11, 12].

Derivations of the different formulas from the different steps in the locks were considered (Figure 2). LS (“left side”) in Figure 2, can represent either the Pacific ocean, the Caribbean Sea, or the volume of water contained in a lock at a lower level. Similarly, RS (“right side”) represents the volume of water contained in a lock at a higher level such as the lakes Miraflores or Gatun. The initial conditions (𝑡=0) for the volume of the lock (𝑉(0)) and the salinity (𝑆(0)) are given by: 𝑆(0)=𝑆0,𝑉(0)=𝑉𝑇𝑉𝐿,(1) where 𝑉𝑇 is the volume capacity of the lock, 𝑉𝐿 is the volume necessary to reach the next level (RS), and 𝑆0 is the initial salinity of the water inside the lock before the water inside and outside are combined. There are four separate steps in an uplockage or in a downlockage process: (1) the gates on the LS lock open, and the water inside the lock combines with the water pushing the ship into the lock; (2) the gates on the LS lock close and the ship is raised to the next water level; the water of both the LS (𝑉LS) and RS (𝑉RS) locks combine and the salinity changes to the new water mix; (3) the gates in the RS lock open, and the water inside the lock combines with the water in the next level (either a lock or a lake); (4) the ship leaves the lock and the water in the lock drains into the lower level (either a lock or the sea); the volume inside the lock is again given by (1).

After the four steps above, the salinity of the water is given by: 𝑆𝑉𝑆(3)=(2)(2)𝑉RS+𝑆RS𝑉RS+𝑉𝑆𝑉(3),(2) where 𝑆(3) is the salinity in steps 3 and 4, since they have the same value. 𝑉𝑆 is the displacement volume of a ship (average) and 𝑉(3) is the volume in step 3. 𝑆(2) y 𝑉(2) are the salinity and volume of step 2, respectively. Equation (2) needs to be applied for each chamber using the known historic salinities as well as the estimated salinity value of the lower level water given by (3): 𝑆𝑁[]=𝑆(0)+𝑁𝑆(4)𝑆(0),(3)𝑆𝑁 provides the salinity value by lock (e.g., Miraflores, Pedro Miguel, or Gatun locks) for any quantity 𝑁 (number of lockages—supposedly only one ship passes through each lockage, so the number of lockages should be the same as the number of ships).

The six locks have different volumes and geometric characteristics so that ships of different drafts can cross the Panama Canal from the Pacific ocean to Gatun Lake. The levels shown in Figure 3 are used with the dimensions shown in Table 1 to calculate the volumes and salinity of the locks when ships cross. In Figure 3 all levels, heights, and depths are referenced to the “Precise Panama Canal Level Reference” (PLD) that matches sea level.

Using the equations of exchange of salinity in the locks, you can set a numerical and differential equations model to define the salinity in Gatun Lake (𝑆GL), taking into account the exchange of water (and salinity) in the upper locks of Pedro Miguel and Gatun, and the water contribution by lakes Gatun (𝑉GL) and the volumetric inflows of Madden (𝑉Madden) with its respective salinity level (𝑆Madden) and the river tributaries (𝑉trib) that flow into these lakes. This relationship is expressed by the following equation:𝑑𝑆GL=𝑉𝑑𝑡GL𝑆GL+𝑉Madden+𝑉trib𝑆Madden+𝑉𝐿3𝑉𝑠𝑆𝐿3𝐸𝑋𝐿3+𝑉𝐿4𝑉𝑠𝑆𝐿4𝐸𝑋𝐿4𝑉GL+𝑉Madden+𝑉trib+𝑉𝐿3𝑉𝑠𝐸𝑋𝐿3𝑉𝑁+𝐿4𝑉𝑠𝐸𝑋𝐿4,𝑑𝑉𝑁GL=𝑉𝑑𝑡Madden+𝑉trib+𝑉𝐿3𝑉𝑠𝐸𝑋𝐿3+𝑉𝐿4𝑉𝑠𝐸𝑋𝐿4Evaporation(𝑡)+Precipitation(𝑡),(4)𝑉𝐿3 and 𝑉𝐿4 are the volumes of Locks L3 and L4, respectively (See Figure 3). 𝑆𝐿3 and 𝑆𝐿4 are the respective salinities of L3 and L4 taking into consideration the measured salinity gradients. 𝐸𝑋L3 and 𝐸𝑋L4 are the exchange ratios for Locks L3 and L4, respectively. Piecewise linear profiles of Evaporation (Evaporation (𝑡)) and precipitation (Precipitation (𝑡)) are added to the calculations of 𝑉GL [1317].

The different equations were implemented in AnyLogic [7, 1820]. AnyLogic provides a set of numerical methods for solving differential equations, algebraic-differential equations, or algebraic equations. AnyLogic chooses the numerical solver automatically at runtime in accordance to the behavior of the system. When solving ordinary differential equations, it starts integration with, fourth-order Runge-Kutta method with fixed steps. Otherwise, AnyLogic plugs in another solver—Newton method. This method changes the integration step to achieve the given accuracy. However, one of the advantages of AnyLogic is its unique hybrid simulation modeling capabilities. We are able to combine the differential equations with a discrete-event system without having to write an interface. This is very important because the diffusion of salinity is a continuous process; however, the lockages produced by the transits of the ships through the canal are essentially a discrete-event process.

3. Discrete-Event Model Operations

In discrete-event simulation, the operation of a system is represented as a chronological sequence of events. Each event occurs at an instant time and marks a change of state in the system. For example, if the Panama Canal is simulated, an event could be a ship entering L1, with the resulting system state. The different steps of the transit can be modeled and the “transit rules” for the respective decisions and the utilization of resources (e.g., locks, lakes) and eventually reaching L6 and leaving the canal. The modeling of the different queues is very important to the implementation of the sequential logic. In addition, data collection and information gathering is essential in order to generate random variables (e.g., Gatun Lake transit time) depending on the Panama Canal current operations.

We conducted interviews to learn about the operations and to obtain the different distributions and times. Figure 4 describes some of the generic operations data obtained from the surveys and interviews with personnel from the Panama Canal.

A discrete-event model was developed in AnyLogic, with the respective animations, Queues, Switches, Java Classes, and the Enterprise (i.e., discrete-event) Library (Figure 5). The Switches were complemented with Java statements to capture the logic of assignment of locks and the schedule of the Panama Canal. A ship from the Atlantic Ocean enters the Panama Canal waters at the port of Cristobal. If the ship is not scheduled to transit that day, it will drop anchor and wait for its scheduled transit time. Otherwise, the vessel will sail toward Gatun Locks with a Panama chief pilot in control of the vessel during its transit. The chief pilot will instruct the ship’s captain as to the speed and direction of the vessel. The chief pilot also will coordinate with the tug operators, line-handlers, and locomotive engineers as to what assistance they need to provide. The chief pilot remains in contact with the Panama Canal Traffic Control Center (TCC) and each lock tower at all times. The captain relays the chief pilot’s instructions to his/her crew members, who perform the proper maneuver [4].

Figure 5 shows the discrete steps of the modeling of a ship going from the Atlantic (ATLANTIC —source using probabilistic functions to generate the ships) to the Pacific (PACIFIC —sink that represent the end of the entity generated by the ATLANTIC process to represent a ship). These discrete process and events are explained as follows. (1)Entity ship with a specific identification (e.g., 101) and its particular features (e.g., 𝑉𝑆) is generated in ATLANTIC (806139.001.graph). (2) The Ship 101 has to wait in a queue (QueueAtlantic 806139.003.graph) for permission to proceed with the transit in the canal (AtlanticEntrance 806139.004.graph—a switch is implemented in Java that communicates and coordinates with other switches to model the decision-making logic of the Panama Canal). Then, when the permission is received, the ship will navigate the strait (StraitAG 806139.005.graph—Strait from the Atlantic to the Gatun Locks) until it arrives to the Gatun Locks where it will have to wait for its respective instructions (GCAQueue 806139.003.graph). Once instructions are received, the ship proceeds using the lane assigned (by the traffic rules programmed in the GCASwitch1 (806139.004.graph) and the GCASwitch2 (806139.004.graph). This logic was complemented with respective Java code as shown in Figure 5 for GCASwitch1). (3) The Canal has two lanes each with its own set of locks and Ship 101 is allowed to enter only one lane in Gatun Locks. (A “Select Output” (806139.006.graph) routes the incoming ship to one of the two lanes depending on the logic and external conditions of the switches. The process of each lock has its own respective operational times and their statistical variations have been modeled using the information gathered from the interviews. The ship enters the first lock of Gatun Locks (LockGCA11 806139.005.graph or LockGCA21 806139.005.graph) depending on the initial decision at the GCAQueue provided by the logic of the switches (GCASwitch1 and GCASwitch2). Next, the ship enters the second lock (LockGCA12 806139.005.graph assuming LockGCA11 was selected initially) or LockGCA22 806139.005.graph (if LockGCA21 was selected). Then, the ship enters the third and final Gatun Lock (LockGCA13 806139.005.graph or LockGCA23 806139.005.graph depending on the initial lane decision at GCAQueue). (4) Finally, at the end of the Gatun Locks a “Combine” (806139.007.graph) which is a synchronization point that lets one ship proceed only after another ship arrives (one can be coming from one lane and another from the other lane). The ship has to wait for the signal to continue at the respective queue (StraitGPQueue 806139.003.graph). When the signal is received, Ship 101 continues navigating Gatun Lake (StraitGP_GatunLake  806139.005.graph) using statistical distributions to model the navigation time and preparing to transit Pedro Miguel Locks. (5) The ship arrives at the Pedro Miguel Locks where it will have to wait for the respective instructions (PMAQueue 806139.003.graph). Then, the ship can proceed using the set of locks assigned (by the traffic rules programmed in the PMASwitch1 (806139.004.graph) and the PMASwitch2 (806139.004.graph) (this logic was complemented with the respective Java code. The ship enters lane one or lane 2 in Pedro Miguel Locks. The ship enters Pedro Miguel Locks (LockPMA1  806139.005.graph or LockPMA2 806139.005.graph—depending on the initial decision at the PMAQueue provided by the logic of the switches GCASwitch1 and GCASwitch2). (6) Finally, at the end of Pedro Miguel Locks, Ship 101 has to again wait for the signal to continue at the respective queue (StraitPMQueue 806139.003.graph) and when the signal is received, it continues navigating Miraflores Lake (StraitPM_MirafloresLake 806139.005.graph) using the corresponding statistical distributions to model the navigation time in Miraflores Lake. (7) The ship arrives at the Miraflores Locks where it will have to wait for the respective instructions (MLAQueue 806139.003.graph). Then, the ship can proceed using the lane assigned (by the traffic rules programmed in the MLASwitch1 (806139.004.graph) and the MLASwitch2 (806139.004.graph). (8)Ship 101 enters the respective lane in Miraflores Locks., The ship enters the first Miraflores Lock (LockMLA11 806139.005.graph or LockMLA21 806139.005.graph—depending in the initial decision at the MLAQueue provided by the logic of the switches MLASwitch1 and MLASwitch2). Next, the ship enters the second lock (LockMLA12 806139.005.graph (if LockMLA11 was selected initially, or LockMLA22 806139.005.graph if LockGCA21 was selected). This second lock is the final Miraflores Lock. (9) Finally, at the end of the Miraflores Locks Ship 101 has to wait for the signal to continue at the respective queue (StraitMPacificQueue 806139.003.graph) and when the signal is received, it continues navigating the strait between Miraflores Locks and Panama Bay (StraitMPacific 806139.005.graph) on its way to the Pacific Ocean. (10) Entity Ship 101 exits at PACIFIC (806139.002.graph) and statistics are generated.

A similar process is required for ships transiting from the Pacific to the Atlantic (where PACIFIC generates ships and ATLANTIC exits the entities). The set of switches and queues simulate the operational logic to allow for the synchronization and coordination needed for these northbound (Pacific/Atlantic) and southbound (Atlantic/Pacific) operations.

4. Hybrid Model

AnyLogic has the capability for hybrid modeling. We combined the system dynamics model with the discrete-event model into one model using the capabilities of the Active Object Class from AnyLogic [7, 19]. The discrete-event model feeds the number of lockages at specific discrete even times and a set of differential and algebraic equations is running in continuous time modeling the salinity diffusion produced by the lockages. The Active Object Class from AnyLogic may have multiple concurrent activities that share object local data and object interface. Activities can be created and destroyed at any moment of the model execution. An activity can be described by a Java function or by a Hybrid statechart. Hybrid statecharts can associate a set of differential and algebraic equations with a simple and/or composite state of a statechart, and the user can also specify a condition over continuously changing variables as a trigger of a transition. The currently active set of equations and triggers is defined by the current simple state and all its containers. Therefore, we can have a discrete-event process and the continuous process running concurrently, sharing information, and influencing their behavior.

5. Animation

The animation of the model was created using Java. The user can watch the ships pass through the locks as well as visualize the salinity changes in the Gatun Lake. The animation is scalable and hierarchical and it gives us an overview of the Panama Canal process and some aggregate indicators such as salinity. A snapshot of the animation is depicted in Figure 6.

6. Validation

The hybrid model was validated using actual transits and salinity field measurements: (i)actual transit data (Years: 2000–2009); (ii)salinity level for the Gatun Lake from field measurements (Years: 2003–2009).

The discrete-event model was developed from interviews with the Panama Canal chief pilots and data obtained for each event. Our validation was performed with actual Panama Canal transit data. In addition, the simulation model was executed by generating 100 sets of independent replications (using different random numbers). We compared the simulation to the system by constructing a confidence interval (CI) as shown in Table 2.

The Salinity level output of the hybrid model was compared with those from the real system (Gatun Lake). The values from the real system were provided by field measurements conducted by the Panama Canal Authority at several points within Gatun Lake. Figure 7 provides the comparison of the hybrid model again the historical data (2003–2009).

In addition, the salinity output of our model is compared against two other models for salinity. The model developed by the Army Corps of Engineers predicts the salinity in the years 2003–2009 to be stable at a value of 0.032 ppt [9]. The Army Corps of Engineers did not consider precipitation, evaporation, hydropower plants, fresh water facilities, and spillway flows. The other model was commissioned by the Panama Canal Authority to Delft Hydraulics (http://www.wldelft.nl/) [12]. Their model does not consider precipitation and evaporation and uses a constant rate of 36 ships/day. The summary of these comparisons is provided in Table 3.

7. Conclusions and Further Research

This research provides a unique example for applying hybrid modeling. Hybrid modeling can benefit organizations with complex operations by providing them with a high prediction capability, which takes into account the internal and external changes taking place in their environment where continuous and discrete variables are present.

The second phase of this project will model the expansion with a third lane and set of locks. Figure 8 indicates the transformation to the current channel. Additions to the current model will have to be performed in order to describe the new traffic patterns involving post-Panamax ships. A lock facility will be located at the Atlantic end of the canal, on the east side of Gatun locks (Figure 8—Point 2). The other facility will be located at the Pacific end of the canal, to the southwest of Miraflores Locks (Figure 8—Point 6) [22].

In addition, the modeling will include new equations to model the salinity diffusion due to the new smart mechanisms being used to build the new locks (water reutilization basins) ensuring the conservation of fresh water. Figure 9 shows the technology to be used for the new set of locks.

Acknowledgment

This work was sponsored by the Secretaría Nacional de Ciencia, Tecnología e Innovación (SENACYT) of the Republic of Panama under Contract/Project no. FID09-079.