*Article* **Combining Optimization and Simulation for Designing a Robust Short-Sea Feeder Network**

#### **Carl Axel Benjamin Medbøen, Magnus Bolstad Holm, Mohamed Kais Msakni \*, Kjetil Fagerholt and Peter Schütz**

Department of Industrial Economics and Technology Management, Norwegian University of Science and Technology, 7491 Trondheim, Norway<EMAIL_ADDRESS>(C.A.B.M<EMAIL_ADDRESS>(M.B.H<EMAIL_ADDRESS>(K.F<EMAIL_ADDRESS>(P.S.) **\*** Correspondence<EMAIL_ADDRESS>
Received: 5 October 2020; Accepted: 16 November 2020; Published: 20 November 2020

**Abstract:** Here we study a short-sea feeder network design problem based on mother and daughter vessels. The main feature of the studied system is performing transshipment of cargo between mother and daughter vessels at appropriate locations at sea. This operation requires synchronization between both types of vessels as they have to meet at the same location at the same time. This paper studies the problem of designing a synchronized feeder network, explicitly accounting for the effect of uncertain travel times caused by harsh weather conditions. We propose an optimization-simulation framework to find robust solutions for the transportation system. The optimization model finds optimal routes that are then evaluated by a discrete-even simulation model to measure their robustness under uncertain weather conditions. This process of optimization simulation is repeated until a satisfactory condition is reached. To find even better solutions, we include different performance-improving strategies by adding robustness during route generation or exploiting flexibility in sailing speed to recover from delays. We apply the solution method to a case based on realistic data from a Norwegian shipping company. The results show that the method finds near-optimal solutions that offer robustness against schedule perturbations due to harsh weather. They also highlight the importance of considering uncertainty when designing a short-sea feeder network with transshipment at sea.

**Keywords:** maritime transportation; liner network design; synchronization; weather uncertainty; optimization simulation

#### **1. Introduction**

Maritime transportation is one of the most efficient (per cargo ton-mile) transportation modes to transport large volumes of cargo over long distances [1,2]. Over short distances, however, competing with road-based transportation is more challenging. Truck-based door-to-door deliveries are often less expensive compared to other shipping solutions and offer frequent and reliable departures [3,4].

The demand for cargo transportation in Norway is expected to grow by 40% (in tonne-kilometers) until 2030 [5]. This predicted increase i is met by the political ambition to shift more goods from road to sea [6]. Still, new and innovative solutions need to be developed to substantially improve the competitiveness of short-sea shipping and support the transition of cargo from road-based to waterborne transportation. The Short Sea Pioneer (SSP) logistics system is a new suggested solution to improve competitiveness, proposing a new transshipment mode for short-sea feeder networks [7]. The proposed system is inspired by the ship-to-ship cargo transfer method currently used to transfer petroleum products and bulk cargo between seagoing vessels. In the SSP system, a mother (large) vessel can be connected to a daughter (small feeder) vessel at a suitable location at sea to transship containerized cargo using a specialized handling cargo system (see [8] for an illustration). One advantage of this system is to reduce the number of port calls, which reduces the operational cost of the shipping system. Indeed, cargo-related port costs can account for up to 30% of the turnover of a smaller short-sea shipping company [9]. Besides, it becomes possible to serve small ports that large vessels cannot visit, for example, due to physical limitations.

The SSP logistics system is composed of one main route sailed by mother vessels that transport cargo between the European continent and large ports located at the Norwegian west coast. Small daughter vessels operate feeder routes serving smaller Norwegian ports. Due to size or location, these ports may only be served by the smaller daughter vessels. Potential candidate transshipment locations will be in sheltered locations (e.g., inside a fjord or inshore) such that harsh weather does not affect the transshipment operation. An example is illustrated in Figure 1.

**Figure 1.** An example of the short sea pioneer (SSP) logistics system [10] composed of one main route and three daughter routes. The main route serves the continental port, Maasvlakte, a large Norwegian port in Bergen, and three ocean hubs. All daughter routes depart from ocean hubs and serve local Norwegian ports.

While potentially having considerable economic benefits, transshipment at sea also raises many technical challenges. Among others, the system requires synchronizing the main route sailed by the mother vessel with all routes sailed by the daughter vessels. Conversely, Weather conditions are known

to impact sailing times and can cause considerable delays [2]. Thus, the synchronization operation can be subject to disruptions that may cause additional waiting times at the transshipment location and can propagate through the system for the subsequent periods. To function correctly, the routes for both mother and daughter vessels need to be robust to potential delays, i.e., they have to account for uncertainty in weather conditions.

Traditionally, the literature on liner shipping network design problems has been focusing on deterministic deep-sea (long-haul) shipping networks. These problems usually determine the optimal set of routes to be served by a heterogeneous fleet of vessels while satisfying demand, transshipment, and frequency requirements (see, e.g., Brouer et al. [11] ). The following publications are recent examples, discussing different variants of the traditional liner network design problem. Meng and Wang [12] include the repositioning of empty containers in the shipping operations when determining the optimal network design. Reinhardt and Pisinger [13] consider more advanced route structures, such as butterfly routes, because they allow for better use of vessel resources. Brouer et al. [11] study one of the largest networks operated by a major liner-shipping company, where cargo can be transshipped several times to take more than one route to be delivered. Karsten et al. [14] extend this problem to include transit time restrictions. Balakrishnan and Karsten [15] incorporate limitations on number of transshipments for the cargo. For the short-sea shipping, Msakni et al. [16] study the impact of different network designs for a local liner shipping company. Fadda et al. [17] address the problem of a roll-on roll-off liner service that operates using a hub-and-spoke network design. Akbar et al. [18] provide an economic analysis of introducing autonomous vessels in short-sea shipping. For a more detailed overview of the literature on liner shipping network design, please also see the surveys by Meng et al. [19], Brouer et al. [1] and Christiansen et al. [20].

Synchronization has received relatively little attention in the literature on maritime transportation. Cargo is usually transshipped in ports, often requiring a sequence of arrivals or specifying a time window. The work of Agarwal and Özlem Ergun [21] is one of the fewest papers to address synchronization directly. The problem considers a combined vessel scheduling and cargo routing problem, where transshipment of cargo is only possible when two routes meet at the same port on the same day. In another work, Andersson et al. [22] study a problem from project shipping where different cargoes may require a synchronized delivery. In land-based transportation, synchronization issues are more common. The reader is referred to Drexl [23] for an overview of the literature on VRPs with multiple synchronization constraints.

The research on uncertainty in maritime service networks mainly distinguishes between two types of uncertainty. The first type is uncertainty in service times, usually port times or sailing times, whereas the second type considers uncertainty in demand. When considering uncertainty in operations, the research focuses on keeping a designed schedule, i.e., satisfy frequency requirements or pickup and delivery time windows. As an example of this work line, the problem examined by Wang and Meng [24] of designing liner vessel routes with uncertainty in port operation times. The uncertainty is related to sailing times as a consequence of making up for the delays. In Song et al. [25], a multi-objective liner shipping service problem with uncertain port times is studied. One of the objectives is to minimize schedule unreliability, which is the probability of the vessels arriving after the scheduled time windows. Conversely, Li et al. [26] study how vessels can recover from delays caused by regular uncertainties and unexpected disruptions. In their work, regular uncertainties may happen both at sea and in port, and their characteristics can be estimated using historical data. When it comes to demand uncertainty, the research usually focuses on designing a maritime transportation system such that demand can be served. For example, Ng and Lin [27] study the problem of fleet deployment under incomplete demand information. Lo et al. [28] present a model for designing a ferry network given uncertain demand. An and Lo [29] study the design of more general transit networks under uncertain demand. However, none of these studies consider uncertainty in service times.

One approach to deal with the uncertainty is to combine simulation and optimization models to provide robust solutions. Fischer et al. [30] use simulation and optimization to evaluate the robustness of tactical fleet deployment plans for roll-on roll-off liner shipping with respect to random disruptions at the operational level. Castilla-Rodríguez et al. [31] study the quay crane scheduling problem in a port terminal and consider the uncertainty from the availability of some delivery vehicles and disruptions in quay crane operations. The authors use simulation-optimization to produce robust quay crane schedules. Layeb et al. [32] develop a simulation-optimization method for multimodal freight transportation systems where the uncertainty is related to demand and travel times. Poeting et al. [33] combine a metaheuristic with a discrete-event simulation to provide robust solutions at the operational level in parcel transshipment terminals.

This paper considers the first type of uncertainty, more specifically uncertain sailing times due to harsh weather conditions. However, in contrast to the papers mentioned above, it studies a short-sea liner network instead of a deep-sea network. A common assumption in papers on deep-sea liner network design with uncertain port times is that vessels can make up for delays by increasing their speed between ports (see, e.g., Wang and Meng [24]). This can be reasonable when legs between ports are long. However, for short-sea networks, the sailing legs are short, and weather conditions may prevent the vessel from reducing or even eliminating a delay.

The problem presented in this paper is an extension to the deterministic Short-sea Liner Network Design Problem with Transshipment at Sea (SLNDP-TS) introduced by Holm et al. [10]. It contributes to the research literature on liner shipping network design in three ways. Firstly, we take into account uncertain sailing times due to harsh weather in the network design problem. Secondly, we consider transshipments at sea, which requires that the synchronization of routes to ensure that mother and daughter vessels are at the same location at the same time. Thirdly, we apply an iterative solution approach combining optimization with discrete event simulation to handle both uncertainties in sailing times and synchronization requirements.

Our solution method is based on the hybrid optimization-simulation method proposed by Acar et al. [34] with the difference of using a discrete-event model instead of an optimization model for the simulations. Unlike commonly used probabilistic models, we use wave height from historical weather data to adequately capture the effect of harsh weather conditions on sailing speed. The simulation model links wave height to vessel speed and evaluates the impact on the synchronization operations. The simulation results are returned to the optimization model to select routes that are less likely to generate delays. Additionally, we develop and test different performance-improving strategies to enhance the solutions' robustness to possible delays. The computational study results show that our approach provides more robust solutions than deterministic ones without a significant increase in costs.

The remainder of this paper is structured as follows. In Section 2, the problem is described in more detail. Section 3 outlines the solution approach. The optimization and simulation models are presented in more detail in Section 4. Results from the computational study are presented in Section 5. We conclude in Section 6.

#### **2. Problem Description**

The problem of designing a short-sea liner network with transshipment at sea has been introduced by Holm et al. [10]. The authors provide a deterministic problem formulation as well as a solution method based on a priori route generation. In this section, we first present the problem before discussing the impact of weather conditions on shipping operations in such a network.

#### *2.1. Network Design Problem*

The shipping company serves each port on a weekly basis. The Norwegian ports are classified according to their size. Small ports are only served by daughter vessels, whereas large ports can be visited by either a mother or a daughter vessel. All mother vessels start their voyage at a European continental main port. Transshipment of cargo between mother and daughter vessels occurs at the so-called ocean hubs, which are suitable locations along the Norwegian coast and offer enough stability for vessels to perform transshipment at sea. For presentation and modeling purposes, each ocean hub is artificially split into a north-going and a south-going ocean hub. A mother vessel serves a north-going ocean hub during its northbound journey and a south-going ocean hub on its southbound journey.

The feeder routes taken by daughter vessels are referred to as daughter routes. Each daughter route is served by one daughter vessel. These routes have a maximum duration of one week to ensure weekly port visits. The route served by the mother vessel is referred to as main route. Since the duration of the main route is typically more than one week, the number of deployed mother vessels is equal to the number of weeks rounded up to the nearest integer, thus ensuring weekly service.

The major activity of the shipping company case is to transport cargo between the main continental port and Norwegian ports. There is also some local demand between Norwegian ports, but can be considered as negligible in this study. The aim is to determine the optimal main and daughter routes and the optimal fleet of mother and daughter vessels to be deployed. This problem is at a tactical level, where the established routes last for typically four to 12 months, and, therefore, the weekly demand is assumed to be known and constant.

Transshipment of cargo between mother and daughter vessels is only possible at ocean hubs. During its northbound journey, mother and daughter vessels meet at a north-going ocean hub where cargo is delivered to the daughter vessel. The same mother and daughter vessels meet again at the corresponding south-going ocean hub during the southbound journey of the mother vessel to transship cargo from the daughter vessel. In the case an ocean hub is the northernmost point of the main route, there is no artificial distinction between north- and south-going ocean hub, and the ocean hub is therefore only visited once.

A solution to the problem consists of a set of main and daughter routes serving the ports and an allocation of vessels to routes (as illustrated in Figure 1). This also includes determining the number of deployed mother and daughter vessels. The mother vessels are considered to have enough capacity to transport all cargo. The size of the daughter vessels can be selected from a given set of available capacities. The problem is separable in the size of the daughter vessel because all daughter vessels must have the same capacity due to the technical requirements of the SSP design. The objective is to minimize the weekly operating costs, including weekly time charter costs for each vessel in the selected fleet, bunker, port, and cargo handling costs.

#### *2.2. Impact of Harsh Weather*

A transshipment in an ocean hub requires that both the mother vessel and the daughter vessel have to be present at the same location at the same time. The meeting point is selected in a sheltered location, such that harsh weather will not affect the transshipment operation. However, harsh weather conditions will affect sailing operations as they may force vessels to slow down, causing delays in the vessels' schedules. We refer to the situation of a mother vessel waiting for a daughter vessel in an ocean hub as a synchronization violation. Synchronization violations delay the waiting vessel, potentially delaying later synchronizations and affecting other vessels in the system.

If a vessel is delayed by too much, it might be unable to complete its route within the maximum allowed duration. This is called a duration violation. Duration violations prevent a weekly port visit frequency because the delay is transferred into the next week. For a logistics system with transshipment at sea to be viable in practice, duration violations must be kept at a minimum.

Figure 2 illustrates synchronization and duration violations. Consider a logistics system consisting of a mother vessel, M, and two daughter vessels, D1 and D2. Daughter vessel D1 visits the north- and south-going ocean hubs, 1n and 1s, respectively. The mother vessel and Daughter vessel D2 meet at ocean hub 2 for picking-up and delivering cargo because ocean hub 2 is the northernmost point. For simplicity, Norwegian ports and the continental main ports are not shown.

**Figure 2.** An example to show duration and synchronization violations due to harsh weather. The mother vessel has to wait for 10 h in ocean hub 1s, causing a synchronization violation and a duration violation of five hours.

In the example, the mother vessel is sailing north from the continental main port. In ocean hub 1n, she meets with daughter vessel D1, and after the transshipment, both vessels depart as scheduled. Neither of the vessels is delayed at this point. The mother vessel continues north to ocean hub 2 without any delay. Daughter vessel D2 is not delayed either, and thus both vessels can synchronize as planned in ocean hub 2. In this example, daughter vessel D2 has 10 h of planned idle time before the arrival of the mother vessel. It could, therefore, be up to 10 h delayed and still synchronize with the arriving mother vessel as scheduled.

After leaving ocean hub 2, the mother vessel starts its southbound journey and sails towards ocean hub 1s to synchronize again with daughter vessel D1. However, harsh weather has caused a delay for daughter vessel D1 after departing from ocean hub 1n. Despite arriving on time in ocean hub 1s, the mother vessel has to wait for daughter vessel D1 and experiences a synchronization violation.

The schedule for daughter vessel D1 has planned idle time of five hours in ocean hub 1s before the arrival of the mother vessel. As the vessel is delayed by as much as 15 h, the mother vessel has to wait for 10 h. As the mother vessel only has five hours of planned idle time in its schedule, the next departure from the continental main port is five hours of delay. The synchronization violation in ocean hub 1s has caused a duration violation of five hours for the mother vessel.

Note that daughter vessel D1 does not incur a duration violation, even though it is 10 h delayed when leaving ocean hub 1s. This is because there is enough planned idle time when sailing between ocean hub to 1s and ocean hub 1n to make up for the delay.

#### **3. Solution Approach**

Our iterative solution method is based on the hybrid optimization simulation framework proposed by Acar et al. [34]. The solution method combines an optimization model (also referred to as the master problem) and a simulation model. The role of the optimization model is to select routes for mother and daughter vessels from a set of a priori generated routes, while the simulation is a discrete event simulation model that evaluates the robustness of solutions (a solution is a combination of routes) with respect to uncertain weather conditions. The simulation model results are used to update the costs of the simulated solutions and then faded back to the optimization model. Alternating between optimization and simulation models is repeated until no new improved solution is found. The proposed solution method is illustrated in Figure 3. According to the classification of Crainic et al. [35], the proposed approach is within the alternate simulation-optimization category.

**Figure 3.** An overview of the optimization-simulation framework. From input data, the routes are generated. The master problem and the simulation model are iteratively solved until no new solution is found.

A static set of routes for the mother and daughter vessels is generated using the provided input data. For each of the daughter vessel types, the master problem then selects the cost-minimizing routes from this set by solving an integer programming model (described in more detail in Section 4.1). Afterwards, the combination of routes chosen by the master problem is simulated using historical weather data to estimate the solution's real-world performance (described in more detail in Section 4.2).

As explained in Section 2.2, a duration violation is caused by harsh weather conditions that lead to a delay long enough to prevent the vessel from completing its round trip within its scheduled duration. Synchronization violations can amplify delays as they may transfer the delay of one vessel to another. A duration violation will automatically lead to an initial disruption when starting the next round trip. The magnitude of a duration violation is defined as the number of hours by which a vessel is late compared to the allowed duration of one round trip.

We add a penalty cost based on the magnitude of duration violations from the simulation to the costs of the selected combination of routes. This additional penalty cost makes solutions (or route combinations) prone to duration violations less attractive in subsequent iterations of solving the master problem. When the master problem is solved again with updated route costs, it might choose a new solution with lower costs. The iterations between the master problem and the simulation model continue until all of the selected solutions have been simulated, and no new improving solutions are found. This feedback approach between the master problem and the simulation model allows generating good solutions based on the trade-off between operational costs and robustness (i.e., penalty costs).

Different strategies to provide robustness against disruptions due to harsh weather and/or synchronization violations can be included in the optimization and/or simulation model (discussed more in Section 4.3). Thus, the framework can also be used to evaluate the potential benefit of introducing performance-improving strategies, for example, permitting speed-ups in case of a delay (see, e.g., Fischer et al. [30], Brouer et al. [36]).

#### **4. Optimization and Simulation Model**

This section provides a more detailed description of the optimization and the simulation models used in the proposed solution method. The optimization model is described in Section 4.1, while the simulation model is presented in Section 4.2. Possible performance-improving strategies that can be evaluated using our solution methods are described in Section 4.3.

#### *4.1. Optimization Model*

The optimization model is based on the approach developed by Holm et al. [10] and consists of a route generation procedure and a master problem. The difference is that, in this study, the optimization model incorporates updated costs from the simulation model. For an efficient route generation, a label-setting algorithm (see, e.g., Irnich [37]) is used to generate a priori the routes for mother and daughter vessels. After the routes are generated, the daughter routes are grouped in subsets according to which main route they can be synchronized with.

The label setting algorithm is used to limit the set of routes introduced to the optimization model. Only non-dominated routes are retained, which means that similar routes composed of the same ports but in a different order and with higher operational costs are eliminated. The route generation assumes that vessels sail at their design sailing speed, i.e., no speed-up is allowed. Here, it should be pointed out that some of the deterministically dominated routes might perform better in our setting with uncertain travel time, e.g., due to a sequence of port visits that allows avoiding harsh weather conditions. However, identifying these routes during the route generation procedure is, in general, too computationally expensive.

Before formulating the optimization model, let us first introduce the following notation:

#### Sets

#### Parameters

Binary decision variables

The deterministic optimization problem can be formulated as follows:

$$\min \sum\_{m \in \mathcal{R}^M} \mathcal{C}\_m^M x\_m + \sum\_{d \in \mathcal{R}^D} \mathcal{C}\_d^D z\_d + \sum\_{s \in \mathcal{S}} \mathcal{C}\_s^S y\_s \tag{1}$$

subject to

$$\sum\_{\mathbf{m}\in\mathbb{Z}^{M}}\mathbf{x}\_{\mathbf{m}}=\mathbf{1},\tag{2}$$

$$\sum\_{\mathbf{m}\in\mathbb{Z}^{M}}\mathbf{x}\_{\mathbf{m}}\odot\mathbf{M}\cdot\sum\_{\mathbf{m}\in\mathbb{Z}^{M}}\mathbf{x}\_{\mathbf{m}}\tag{3}$$

$$\begin{aligned} \sum\_{d \in \mathcal{R}\_p^D} z\_d - M^{\text{OH}} \sum\_{m \in \mathcal{R}\_p^M} \mathbf{x}\_m &\le 0, & \mathbf{p} \in \mathcal{P}^{\text{OH}},\\ \sum\_{d \in \mathcal{D}\_p} z\_d - \sum\_{m} z\_d &\ge M^{\text{OH}}(\mathbf{x}\_m - 1), & \mathbf{p} \in \mathcal{P}^{\text{OH}}, m \in \mathcal{R}\_p^M, \end{aligned} \tag{3}$$

$$\begin{aligned} &d \in \mathbb{R}\_{\text{pw}}^{D} \\ &\sum\_{d \in \mathbb{R}\_{\text{pw}}^{D}} z\_d - x\_{\text{m}} \ge 0, \end{aligned} \tag{5}$$

$$\sum\_{\mathbf{m}\in\mathcal{R}\_{p}^{M}} \mathbf{x}\_{\mathbf{m}} + \sum\_{d\in\mathcal{R}\_{p}^{D}} \mathbf{z}\_{d} = \mathbf{1}, \tag{6}$$

$$\sum\_{d \in \mathbb{R}\_p^D} z\_d = 1,\tag{7}$$

$$\sum\_{m \in \mathcal{R}^M} S\_{ms}^M \mathbf{x}\_m + \sum\_{d \in \mathcal{R}^D} S\_{ds}^D \mathbf{z}\_d - \mathbb{S}\_s + \varepsilon \le M^S y\_{s\prime} \tag{8}$$

$$\mathbf{x}\_{\mathfrak{m}} \in \{0, 1\}, \tag{9}$$

$$\mathfrak{m} \in \{0, 1\}, \tag{10}$$

$$y\_s \in \{0, 1\}, \tag{10}$$

$$z\_d \in \{0, 1\}, \tag{1}$$

The objective function (1) minimizes the total weekly cost of the shipping system. The first and second terms of (1) are related to the costs of mother and daughter vessels deployed in the system. The last term represents the penalty cost of a solution.

Equation (2) forces the optimal solution to select only one main route. Constraints (3) to (5) are the synchronization constraints, ensuring that selected daughter routes can synchronize with the main route at an ocean hub visited by both routes. The reader is referred to Holm et al. [10] for a more detailed description of how synchronization is taken care of through these constraints. Further, Equation (6) ensure that a main port is visited by either a mother or a daughter vessel. Equations (7) ensure that each small Norwegian port is served by a daughter route. Constraints (8) set the indicator variable *ys*. If the master problem selects a route combination that constitutes a previously simulated solution, *ys* = 1 and the penalty cost is added to the objective. Lastly, Constraints (9)–(11) restrict the variables to take binary values.

#### *4.2. Simulation Model*

The optimal solution of the master problem assumes that the shipping system runs under perfect weather conditions, and both mother and daughter vessels sail at design speed. The simulation model evaluates the proposed solution under real-world operating conditions. Harsh weather conditions may force vessels to slow down, affecting the sailing times. In such a situation, synchronization and duration violations can occur.

The simulation model adjusts the sailing speed based on wave height. To this end, we divide a route into legs connecting two consecutive ports. For a given leg, waypoints are defined and mark a change in the vessel's travel direction. The segment between two successive waypoints is defined as a sub-leg. The simulation model updates the sailing speed of a vessel at specific points of a leg. Such a point is called an observation point and is defined (i) at a departure port and (ii) after sailing a certain distance, called step distance.

Figure 4 shows an example of a route between two ports. There are three sub-legs defined by the ports and waypoints. The observation points are equally separated by the step distance, which is set to be five nautical miles in this example. The background of Figure 4 shows the gridded weather data points.

**Figure 4.** Illustration of a leg. The vessel changes direction at the waypoints and the speed is updated at every observation point.

To determine the speed of a vessel, we use historical weather data composed of a grid of data points that contain significant wave height. The parameter 'significant wave height' is often used to describe the weather in maritime navigation and is defined as the average height of the highest one-third waves [38]. At each observation point, we extract the significant wave height from the closest data point.

Estimating the vessel speed in open water requires detailed information about the vessel, e.g., regarding hull design and loading conditions, to carry out accurate hydrodynamic calculations. Such information is often not available at an early stage of designing feeder networks. Therefore, we use a simplified approximation of the hydrodynamic relationship between wave height and vessel speed in our simulation model. We assume that vessels are often designed to be operated at a particular propulsion power for which the engines run efficiently, and they sail with constant design power that is maintained throughout the route. Since the vessel speed is a function of propulsion power and total resistance, and given that the propulsion power is assumed to be constant, higher waves and winds due to harsh weather cause an increase in total resistance, which results in a decrease of vessel speed. In our model, resistance from waves is calculated according to the STAwave-1 method [39].

The flow chart in Figure 5 shows the overall simulation process. The input to the simulation model is one (or more) feasible solution(s) from the master problem. In Step 1, the simulation model chooses a solution from the master problem. The week to simulate is chosen in Step 2. In Steps 3–5, all routes in the solution are simulated. Once all routes in the solution are simulated for the given week, synchronization and duration violations are calculated in Step 6. Note that these violations are system-specific, i.e., dependent on the set of routes in the given solution and not only on each route individually. If any violations occur, penalty costs are added to the solution in Step 7. Step 8 checks whether the solution has been simulated for all weeks. Finally, Step 9 ensures that all solutions provided by the optimization model are simulated.

**Figure 5.** Flow chart of the simulation model. The solutions provided by the master problem are simulated week-by-week over a long period using historical data. When a violation is detected, a penalty cost is added to the operational costs of the solution.

An important element of the simulation model is to penalize duration and synchronization violations. The total penalty cost for a simulated solution is the accumulated penalty costs for all vessels for each week of the simulation period. Hence, route compositions that often incur duration or synchronization violations are receiving a high penalty cost and are less likely to be chosen by the master problem in subsequent iterations.

#### *4.3. Performance-Improving Strategies*

Penalizing synchronization and duration violations helps select more robust routes for both mother and daughter vessels, but does not necessarily reflect real-world operations of a feeder network service. In our optimization-simulation model, duration violations propagate from one week to another. In contrast, real-world vessel operators would try to either avoid a delay (e.g., by using more robust routes) or recover from the delay (e.g., by speeding up) before it causes any violations. To further improve the quality of the solutions, we propose two types of performance-improving strategies that mimic real-world operations and are integrated into the solution approach: enhancing the robustness of the routes and recovery from delays.

#### 4.3.1. Enhancing Route Robustness

The first type of strategy tries to enhance the routes' robustness during the route generation procedure, i.e., before starting the optimization-simulation process. Additional slack (or idle time) is added to some or all of its legs when generating routes. By doing so, a vessel might still be able to synchronize in an ocean hub, even if it has to slow down due to harsh weather. Conversely, adding too much slack can cause excessive waiting times and, as such, increase the operational cost for the vessels. Therefore, determining the amount and place of additional slack while ensuring a sufficient robustness level is an important trade-off to assess. This strategy is comparable to the current practice of airlines where delayed routes can recover at night when no flights are scheduled and the network resets [36].

Similarly, the robustness of the routes can be enhanced by considering seasonality in the schedules. Due to the difference in weather conditions between summer and winter, it might be beneficial to operate dedicated routes for each season. However, the system needs to operate the same fleet of mother and daughter vessels during the whole year, regardless of the season.

#### 4.3.2. Recovery from Delays

The second type of strategy tries to exploit operational flexibility to recover from delays. If weather permits, a vessel can increase its speed to reduce or even eliminate a delay. Fischer et al. [30], for example, show that speed-up can be used as a recovery action when considering disruptions in planning fleet deployment in roll-on roll-off liner shipping. Though speeding-up allows recovering from delays in some cases, it requires increased power output, which in turn will increase fuel consumption and thus bunker costs. The ability to speed-up is limited by the maximum propulsion power of the vessels.

#### **5. Computational Study**

The route generation procedure and the simulation model have been implemented in MATLAB R2016b. The master problem is formulated and solved using FICO Xpress 7.9. All computations have been carried out on a PC with a 3.4 GHz Intel Core i7 processor and 32 GB RAM.

#### *5.1. The Short Sea Pioneer Logistics System*

The selection of ports is based on the ports currently served by the case shipping company NorthSea Container Line (NCL), which operates a line for container transportation between the Norwegian west coast and the European continent port located at Maasvlakte in the Netherlands. The candidate hub locations are selected close to major local ports and inshore to have stable weather conditions for the cargo transshipment. The port locations of this case study are illustrated in Figure 6.

The weekly number of containers transported to and from the European continent is based on the number of containers imported and exported in ports along the Norwegian west coast during the first quarter 2016 data [40]. The numbers have been corrected according to NCL's market share to approximate realistic transportation demand. We refer to Holm et al. [10] for detailed input data.

The possible capacities for the daughter vessels are 100, 200, and 300 twenty-foot-equivalent units (TEU). The capacity of the mother vessels is 2500 TEU, which is sufficient to transport the cargo of the whole system. The design speed for all vessels is set to 12 knots.

The relationship between speed and wave height for the mother and three candidate daughter vessel types considered in this experiment is shown in Figure 7. We can see from Figure 7a that at design power, the sailing speed significantly decreases as waves become higher, especially for the daughter vessels. Figure 7b illustrates the same relationship, but for speed-up power that can be used to recover from delays. The power output for speed-up is set to be 125% of design power.

The legs between ports required to simulate sailing speed have been discretized according to container vessel AIS data available from the Norwegian Coastal Administration. Historical wave height data from both the Norwegian Meteorological Institute (MET.no) and the European Centre for Medium-Range Weather Forecasts (ECMWF) is used in the simulations. These data sets provide weather data with a time increment of six hours. Figure 8 shows a visualization of the data grid for both data sources. The data from ECMWF has a resolution of 0.75◦ × 0.75◦, which is sufficient for the zone between Europe and Norway but not detailed enough for vessels sailing the Norwegian coastline. The data from Norwegian Meteorological Institute offers a better resolution, i.e., 0.3◦ × 0.15◦ and is used for latitudes higher than 57.5◦ N. Please note that no wave height data is available inside fjords or very close to land. If a vessel is within these areas, it is assumed to sail at design speed, i.e., 12 knots.

In our simulations, we use data of two successive years (2000 and 2001). To make the simulation by weeks easier, we assume that the length of each month is exactly four weeks, resulting in a total of 96 weeks to simulate. The winter season is defined as lasting from October to March, and the summer season lasts from April to September.

**Figure 6.** A map showing ports considered in the case study. The possible ocean hubs are represented with yellow squares.

**Figure 7.** Impact of significant wave heights on the sailing speed for the fleet of vessels considered in the experimentation.

(**a**) European Centre for Medium-Range Weather Forecasts (ECMWF) data used south of a latitude of 57.5◦N

(**b**) Data from MET.no used north of 57.5◦N

**Figure 8.** Data grid with time series observation on each cell used for the weather simulation.

#### *5.2. Runtime Considerations*

The runtime of the optimization-simulation framework mainly depends on the number of iterations needed for the framework to converge to a solution. To speed up this process, we can reduce the number of times to run the optimization model by returning multiple feasible solutions for simulation instead of only the optimal solution. This has the advantage of possibly simulating the solution that performs best under uncertainty at early iterations. On the other hand, extracting too many solutions implies simulating solutions that may never be needed. Thus, it would be interesting to study the trade-off between the number of solutions extracted in each iteration and the number of iterations of the optimization-simulation framework.

We extract the *n* best solutions found while solving the master problem. The results for the convergence speed of the solution method for different values of *n* are presented in Table 1. The first column indicates the maximum number of extracted solutions *n* in each iteration. The column Itr. needed states how many times the master problem has been solved, i.e., the number of iterations. The third column, Route Gen. Time, is the time needed to generate the routes. Note that this time is constant, as the routes are only generated once. The total time spent by Xpress to solve the master problems is given in column Opt. Time, whereas the total time used by MATLAB to simulate the solutions is given in column Sim. Time. Column Total Time shows the total runtime spent by the solution approach. The last column, Solutions Simulated, contains the total number of simulated solutions.

As can be seen from Table 1, the number of iterations of solving the master problem decreases with higher values of *n*. In contrast, the time needed to simulate all solutions increases by *n*. For our case, we select *n* = 20 for which the best total time is obtained. Finally, we point out that the approach is not sensitive to *n* when considering the solution quality. The same best solution is found for all different values of *n*.

**Table 1.** Runtime (in seconds) for the solution approach with different values of *n*.

#### *5.3. Determining Global Estimated Optimal Routes*

In this section, we analyze the route networks obtained using the solution method. We first describe the different performance-improving strategies, before discussing the solutions produced by our method.

#### 5.3.1. Performance-Improving Strategies

We combine the optimization-simulation solution method with different performance-improving strategies to impose additional robustness and exploit operational flexibility. The different strategies are listed in Table 2 and are explained in more detail below.

We test three different performance-enhancing strategies for enhancing the robustness of the routes and one strategy for recovery from delays. The Slack 10% and the Realistic strategies are implemented as part of the route generation procedure. The Slack 10% strategy is characterized by calculating the sailing time along a leg (based on design speed) and then adding 10% additional sailing to immunize the route against delays. While adding robustness to a route, this strategy may result in considerable amounts of unnecessary idle time in the schedule. The Realistic strategy tries to account for the fact that particular sailing legs are more likely to be subject to delays than others. Sailing speeds along the different legs are first simulated using weather data from 1998 to 2000. To further improve robustness, the resulting average sailing time on each leg is then increased by an additional 5%.

With the Seasonal strategy, the routes sailed by the vessels are allowed to change for each season. We distinguish here between the summer season and winter season. The seasonal routes are obtained by solving the simulation-optimization framework for each season with a set of routes based on design speed. Note that the same fleet has to be used in all seasons.

The Speed-up strategy is part of the simulation model. This strategy allows delayed vessels to increase their speed by increasing their power output. Mother and daughter vessels speed-up independently according to their routes, but the conditions for when they do so are slightly different. A mother vessel will try to speed-up once it is more than one hour delayed compared to its generated route. Increased speed will be maintained until the vessel is no longer delayed. Daughter vessels usually have a certain amount of idle time included in their schedules, e.g., waiting for the mother vessel to arrive in an ocean hub. The daughter vessel starts to speed-up as soon as the remaining idle time falls below a given threshold (for the case study, this is set to one hour) and will continue to sail

at increased speed until it has up to two hours idle time available. However, a daughter vessel will never speed-up to get more than its original amount of idle time.

The Combined strategy applies realistic, speed-up, and seasonal strategies simultaneously. The realistic sailing times used in the combined strategy are calculated separately for each season during the preliminary simulation of sailing time. This implies that more absolute slack will be added in the winter season than in the summer season.

#### 5.3.2. Best Routes with and without Weather Uncertainty

We first solve the problem deterministically and use its solution as a benchmark when comparing the solutions considering weather uncertainty. For the deterministic case, the solution using daughter vessels with a capacity of 300 TEU has the lowest total weekly cost with an objective function value of 457,800 USD. Two mother vessels and two daughter vessels are deployed in this optimal solution, as shown in Figure 9. In particular, the main route visits two ocean hubs at Tananger and Haugesund, as well as the main port at Tananger. The hub Haugesund is the northernmost port served by the main route and is therefore visited only once. The mother vessel is more expensive in terms of bunker cost and port visits, which causes the model to select a short main route.

**Figure 9.** The optimal routes for the Short Sea Pioneer (SSP) deterministic version. The mother route visits two hub ports that are connected to two daughter routes.

We then solve the network design problem using our optimization-simulation method, applying each of the performance-improving strategies separately. In all best solutions, two mother vessels and two daughter vessels are deployed. However, the solution without the performance-improving strategy deploys smaller daughter vessels than other solutions. Table 3 summarizes the results for the different performance-improving strategies. The cost of the best solutions is compared to the cost of the best deterministic solution, where 100% representing the cost of the deterministic solution. In this comparison, the penalty costs added by the simulation

model are excluded as they are fictional and only used to identify low-performing solutions. However, the increased costs due to speeding up are included. This comparison facilitates evaluating the solution quality since the deterministic solution is obtained under ideal weather conditions and, thus, can be considered as a lower bound. Table 3 also shows important performance characteristics for the different solutions, such as the total number and accumulated time (in hours) of both duration and synchronization violations. It also reports the total amount of planned idle time included in the selected routes and how much time this idle time is used to recover from incurred delays. The last column indicates the daughter vessel size of the solutions.

**Table 3.** Performance characteristics for the different solutions.

All obtained solutions result in a slight increase in operational costs compared to the deterministic solution, where the Combined strategy gives the lowest operational cost with an increase of only 1.2%. Only one duration violation occurs with a negligible magnitude. The number of synchronization violations and the related magnitude is also significantly lower than for other strategies. The utilization of idle time is also low, implying that the vessels usually keep the estimated arrival times set up in the route generation procedure. As such, the solution obtained when using the combined performance-improving strategy is quite robust without increasing the cost level by much. Note that the solution resulting from not using a performance-improving strategy does not incur any duration violation but is exposed to many synchronization violations. The solution includes a lot of idle time, a large share of which needs to be used for sailing to recover from delays.

Figure 10 shows the optimal seasonal routes for the Combined strategy. The solid red line is the main route that continues further south to Maasvlakte port. The blue and pink dashed lines are daughter routes 1 and 2, respectively. The corresponding ocean hubs are marked as squares.

As seen in Figure 10, the main route extends further north in the winter than during the summer season. Compared to a daughter vessel, a mother vessel is more expensive in terms of bunker and port costs. By having a mother vessel sailing further north, the daughter vessels can reduce their total sailing distance, resulting in planned idle time. This is beneficial during the winter season since the weather conditions are rougher, and a greater buffer against delays is needed to avoid duration violations. Conversely, during the summer season, the weather conditions are better, resulting in a shorter main route with corresponding longer daughter routes to reduce operational cost.

#### 5.3.3. Analyzing the Best Deterministic Solution

The best deterministic solution is slightly cheaper than all solutions found by the solution triggered feedback approach. However, these cost savings come at the expense of a lack of robustness for harsh weather conditions. Table 4 provides the performance characteristics for simulating the optimal deterministic solution under weather conditions without any performance-improving strategy.

(**a**) Routes for the summer season (**b**) Routes for the winter season.

**Figure 10.** The difference between the summer and winter solutions when using the combined strategy. The main route sails further north in the winter season due to rough weather conditions.

In particular, the high number of duration violations renders the best deterministic solution infeasible in practice. While almost all of the planned idle time is used for sailing to mitigate delays, the system is still unable to recover. The poor performance of the optimal deterministic solution clearly shows the importance of taking into account weather uncertainty. All solutions found the using optimization-simulation approach (see Table 3) perform significantly better, while being less than 3% more expensive. The solution from the Combined strategy for example, reduces the number of duration violations from 277 to 1 and the duration magnitude from 45,429 h to 0.7 h, while the costs only increase by 1.2%.

Pushing the analysis further, we simulate the best deterministic solution with the Speed-up strategy. This combination results in a cost increase of 2.2%, which is due to the higher fuel consumption from the required speed-ups. However, the number of duration violations is reduced from 277 to 28, and the number of synchronization violation is reduced from 194 to 181. Even though these are significant improvements, the solution still performs much worse than the solutions found by the optimization-simulation method.

In another analysis, we consider another criterion of a well-performing solution. It consists of measuring the stability of arrival times, which translates directly into predicting future arrivals and maintaining a given schedule. For the Combined-strategy solution, the arrival times for the summer and winter routes are presented in Figures 11 and 12. The mother vessels are scheduled to depart from the main continental port at time zero. The arrival times for all vessels refer to this departure time.

**Figure 11.** A box and whisker plot of the arrival times for the summer season of the Combinedstrategy solution.

**Figure 12.** A box and whisker plot of the arrival times for the winter season of the Combinedstrategy solution.

The small vertical line represents the median arrival time. The width of each box corresponds to the interquartile range (*IQR*) defined by the first and third quartile, and accounts for 50% of all port arrivals in the simulation period. The left and right whiskers capture all arrival times within 1.5 · *IQR*. Arrival times outside this range are considered outliers and represented with circles.

For both seasons, arrival times are mostly within relatively small intervals, particularly for the main route and daughter route 1. The longer daughter route 2 is slightly more prone to variations in arrival time, especially during the winter season, although this route is considerably shortened during this season (see Figure 10). Still, these varying arrival times do not affect synchronization with the mother vessel at ocean hubs (except Bergen) because daughter vessel 2 always arrives first anyway (with one exception).

Figure 13 presents the arrival times for the deterministic solution. Here, we observe much wider intervals for different arrival times. For example, about 50% of all arrivals at the main continental port in Maasvlakte happen approximately 48 h around the median arrival time. This median arrival time of the mother vessel is also about 100 h after the end of the planned two-week route due to a large number of duration violations. Moreover, Figure 13 shows that the deterministic solution is very sensitive to weather uncertainty, causing arrival times to vary a lot. Therefore, the deterministic solution is not applicable in a real-world setting and should be discarded in favor of the solution from the Combined-strategy. These results highlight once more the need for considering the weather uncertainty in designing the short-sea feeder network considered in this paper.

**Figure 13.** A box and whisker plot of the arrival times of deterministic solution.

#### **6. Concluding Remarks**

In this paper, we study the problem of designing a short-sea liner network with transshipment at sea under uncertain weather conditions and propose an iterative solution method that combines optimization and simulation. In the deterministic optimization step, we select routes for both mother and daughter vessels that minimize the cost of operating the logistics system. These routes are then evaluated using a discrete event simulation model. Solutions that are not performing well in the simulations are penalized in the objective function in subsequent iterations of the optimization model to reduce their attractiveness. The steps of selecting routes in the optimization model and simulating them continue until no new (not simulated) solution is selected. The computational study performed on a case based on real-world data shows that our solution method provides solutions for the system with a good trade-off between operational cost and robustness.

The results show that weather uncertainty can severely impact the synchronization of the routes and should be taken into account in the design phase of the logistics system. The optimization-simulation approach, especially when using different performance-improving strategies, finds robust solutions at only a small operational cost increase. The solution method has been applied to a short-sea feeder network design problem in this paper, but it should be straightforward to apply to other related problems.

**Author Contributions:** Conceptualization, K.F. and P.S.; methodology, C.A.B.M., M.B.H., K.F. and P.S.; software, C.A.B.M. and M.B.H.; validation, C.A.B.M., M.B.H., K.F. and P.S.; formal analysis, C.A.B.M. and M.B.H.; investigation, C.A.B.M. and M.B.H.; writing—original draft preparation, C.A.B.M. and M.B.H.; writing—review and editing, M.K.M., K.F. and P.S.; visualization, C.A.B.M. and M.B.H.; supervision, K.F. and P.S.; project administration, K.F. and P.S. All authors have read and agreed to the published version of the manuscript. Authors contributed equally to this work.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).


**3.3 Soils**

According to the results we found higher values of PAHs in crops than soil. For the case of crops, the individual PAHs were slightly lower in wet than dry season; in the rainy season the irregular intensity and duration of rainfall cleans the atmospheric of suspended particles. In figure 4 the results show greater values in apple from Tlahuac and cactus stem from Milpa Alta in the wet season. In general, the compounds > 3 aromatic rings were

Presence of Polycyclic Aromatic Hydrocarbons

(PAHs) in Semi-Rural Environment in Mexico City 225

Fig. 5. Values of PAHs in soils from semi-rural terrains from Mexico City.

volatile compound.

as a source of contamination.

Nap not found in our analysis we supposed a loss in handling of sample due that is high

The most prominent source of PAHs in the urban environment is the incomplete combustion of biomass (such as vegetation) and fossil fuels (petroleum). Vehicular traffic (mainly diesel-powered) is considered to be the most significant contributor to the atmospheric PAH load within urban areas (Marr et al., 2004). Atmospheric deposition is the most common source of pollution in soil and it is expected that most combustion derived PAHs will be restricted to the top layer of the soil. Urban areas generally have high traffic density which results in heavy contamination of surface soils (Agarwal, 2009). For our study areas there are several asphalt ways that communicate with population areas, it considering

According Cram et al (2008), the high concentration of these contaminants in atmospheric is concentrate in South of Mexico City, mainly due to the direction of wind from North to South. As appreciate in figure 1 Tlahuac and Milpa are located in Southeast. These soils are considering as conservation, recharge of aquifer and recreational areas further gives other environmental services as retention, regulation and alleviates the deposition of contaminants derived of fossil fuels combustion that origins heavy metals and ethers, hydrocarbons, etc.

present in crops and soil in both season such as BaA, Cry, Fla, Pyr, DaA, BaP and BkF are dominant in apple while cactus stem. According García-Alonso et al. (2003) the presence of Fla and Bghi in the two areas may indicate a common vehicular emission source.

There are many factors influencing the distribution of different PAHs in the environment, such as physicochemical properties of PAHs, physicochemical properties of soil, sources of emission of PAHs, and photochemical degradation of atmospheric PAHs. Two to three ring PAHs are subject to atmospheric transport to remote areas and are considered ''multi-hop'' chemicals, while higher ring PAHs are associated with particles and undergo ''single hop'' transport behavior, and higher ring PAHs are prone to rapid deposition and retention close to source regions. So, PAHs may become fractionated from source regions to remote regions during atmospheric transport (Wang et al., 2010).

The concentrations of BaP in the soils is an indication of both pyrogenic and petrogenic sources of PAH pollution on the environment (Essumang et al., 2010) for both seasons. BaP is considered within permissible limits (2 mg/kg) for agricultural and residential land use according to Mexican regulation NOM 138-SEMARNAT/SS-2003. However, the levels of BaP surpassed the limit established by the Danish Environmental Protection Agency (0.1 mg/kg) and Canadian Council of Ministers of the Environment (0.26 mg/kg; Essumang et al., 2010). The BaP concentrations are considered to be a risk to human health due to its potential exposure to carcinogenic PAHs for those living in those areas. According other international regulations our soils has severe problems of contamination, where some prevent measurements have applied such as strict control of emissions of industries, vehicular park and matter combustion.

Within polycyclic aromatic hydrocarbons considered in Mexican regulation, the concentrations of DaA, BaA, BbF, BkF and Ind were lower than the permissible limit for agricultural land use according NOM-138-SEMARNAT/SS-2003 (2 mg/kg). These compounds do not represent a risk for human health in comparison to BaP.

In figure 5 we appreciated the distribution of PAHs in dry and wet season found it a relative high concentration of chrysene in wet season in both sites. In general dominant PAHs for dry season were BaA, Cry, BaP, Fla, Ind and BkF with range of 0.5 to 3.5 mg/kg; Milpa Alta has slightly high BaA, BaP and Fla and Tlahuac were Cry and DaA. In wet season dominant PAHs were Cry and BaP mainly, with range of 1 to 9.5 mg/kg; Tlahuac has high concentration of Cry and BaA and Milpa Alta were Fla, BaP, Bghi, BbF, BkF, Ind and DaA. The variability of individual PAHs is key to recognize the possible source of contamination in the semi-rural area further the environmental condition defines distribution and concentration of contaminants in the soil and other media.

The presence of 4-6 aromatic rings compounds in soils in Mexico City was similar to results from a survey in soil described by Wang et al. (2010) in urban soils from North China. In semirural soil of Mexico the compounds predominant were Cry, BaP and BaA different at found it in North of China. For dry season we found similar percent of four (50 % approximately); five (25 % approximately) and six (10 %) ring aromatic compounds in Tlahuac and Milpa Alta (Figure 6). In wet season, there are significant changes over aromatic compounds, in Tlahuac has 90 % of four ring aromatic compounds while others compounds diminished to 8 % (5 rings compounds) and 2 % (6 ring compounds). So, Milpa Alta shows similar percent of aromatic compounds than dry season.

present in crops and soil in both season such as BaA, Cry, Fla, Pyr, DaA, BaP and BkF are dominant in apple while cactus stem. According García-Alonso et al. (2003) the presence of

There are many factors influencing the distribution of different PAHs in the environment, such as physicochemical properties of PAHs, physicochemical properties of soil, sources of emission of PAHs, and photochemical degradation of atmospheric PAHs. Two to three ring PAHs are subject to atmospheric transport to remote areas and are considered ''multi-hop'' chemicals, while higher ring PAHs are associated with particles and undergo ''single hop'' transport behavior, and higher ring PAHs are prone to rapid deposition and retention close to source regions. So, PAHs may become fractionated from source regions to remote regions

The concentrations of BaP in the soils is an indication of both pyrogenic and petrogenic sources of PAH pollution on the environment (Essumang et al., 2010) for both seasons. BaP is considered within permissible limits (2 mg/kg) for agricultural and residential land use according to Mexican regulation NOM 138-SEMARNAT/SS-2003. However, the levels of BaP surpassed the limit established by the Danish Environmental Protection Agency (0.1 mg/kg) and Canadian Council of Ministers of the Environment (0.26 mg/kg; Essumang et al., 2010). The BaP concentrations are considered to be a risk to human health due to its potential exposure to carcinogenic PAHs for those living in those areas. According other international regulations our soils has severe problems of contamination, where some prevent measurements have applied such as strict control of emissions of industries,

Within polycyclic aromatic hydrocarbons considered in Mexican regulation, the concentrations of DaA, BaA, BbF, BkF and Ind were lower than the permissible limit for agricultural land use according NOM-138-SEMARNAT/SS-2003 (2 mg/kg). These

In figure 5 we appreciated the distribution of PAHs in dry and wet season found it a relative high concentration of chrysene in wet season in both sites. In general dominant PAHs for dry season were BaA, Cry, BaP, Fla, Ind and BkF with range of 0.5 to 3.5 mg/kg; Milpa Alta has slightly high BaA, BaP and Fla and Tlahuac were Cry and DaA. In wet season dominant PAHs were Cry and BaP mainly, with range of 1 to 9.5 mg/kg; Tlahuac has high concentration of Cry and BaA and Milpa Alta were Fla, BaP, Bghi, BbF, BkF, Ind and DaA. The variability of individual PAHs is key to recognize the possible source of contamination in the semi-rural area further the environmental condition defines distribution and

The presence of 4-6 aromatic rings compounds in soils in Mexico City was similar to results from a survey in soil described by Wang et al. (2010) in urban soils from North China. In semirural soil of Mexico the compounds predominant were Cry, BaP and BaA different at found it in North of China. For dry season we found similar percent of four (50 % approximately); five (25 % approximately) and six (10 %) ring aromatic compounds in Tlahuac and Milpa Alta (Figure 6). In wet season, there are significant changes over aromatic compounds, in Tlahuac has 90 % of four ring aromatic compounds while others compounds diminished to 8 % (5 rings compounds) and 2 % (6 ring compounds). So, Milpa

compounds do not represent a risk for human health in comparison to BaP.

concentration of contaminants in the soil and other media.

Alta shows similar percent of aromatic compounds than dry season.

Fla and Bghi in the two areas may indicate a common vehicular emission source.

during atmospheric transport (Wang et al., 2010).

vehicular park and matter combustion.

Fig. 5. Values of PAHs in soils from semi-rural terrains from Mexico City.

Nap not found in our analysis we supposed a loss in handling of sample due that is high volatile compound.

The most prominent source of PAHs in the urban environment is the incomplete combustion of biomass (such as vegetation) and fossil fuels (petroleum). Vehicular traffic (mainly diesel-powered) is considered to be the most significant contributor to the atmospheric PAH load within urban areas (Marr et al., 2004). Atmospheric deposition is the most common source of pollution in soil and it is expected that most combustion derived PAHs will be restricted to the top layer of the soil. Urban areas generally have high traffic density which results in heavy contamination of surface soils (Agarwal, 2009). For our study areas there are several asphalt ways that communicate with population areas, it considering as a source of contamination.

According Cram et al (2008), the high concentration of these contaminants in atmospheric is concentrate in South of Mexico City, mainly due to the direction of wind from North to South. As appreciate in figure 1 Tlahuac and Milpa are located in Southeast. These soils are considering as conservation, recharge of aquifer and recreational areas further gives other environmental services as retention, regulation and alleviates the deposition of contaminants derived of fossil fuels combustion that origins heavy metals and ethers, hydrocarbons, etc.

Presence of Polycyclic Aromatic Hydrocarbons

and recharge of aquifers with low density of population.

Compounds Soils (µg/kg)

Tláhuac Milpa Alta Fangcun Rural

study, 2009; 2. Chen et al., 2005; 3. Agarwal, 2009; 4. Nam et al, 2003.

employ for irrigation for crops or human and animal consumption

Ingestion of particles or food on contaminate soils

Inhalation of smallest particles

(Sabroso & Pastor, 2004).

**3.4 Sources of PAHs in semirural sites**

The importance of contamination of soil is due serious risks on population health due:

Direct contact with skin in workers of field or direct markets of commercialization

Further the contamination of water resource in form superficial or underground that are

Government hopes that in medium time the contaminants reduce significantly in favour the environment and health population of Mexico by programs of improvement of quality air.

We used the Flu/(Flu+Pyr), Ant/(Ant+Phe), Ind/(Ind+Bghi), BaP/(BaP+Cry), BkF/Bghi and BaP/Bghi ratios and the majority of samples fell into the section identifying pyrogenic sources (fossil fuel, grass and garbage combustion) (Table 3). This is logical considering the

(PAHs) in Semi-Rural Environment in Mexico City 227

industrialisation (Marr et al, 2004; Maliszewska-Kordybach et al., 2009). The spatial distribution of soils contaminated with PAHs reflected that rural/forested/ recreation areas decreasing along the South to North of Mexico City. Due in North there high concentration of industries and urban areas while in South of City the soils are considered conservation

Nap --- --- 131 --- 78.5 Acy --- --- 50 120 317 41.5

Ace 60 180 9 198 298 33.7 36.0 Flu --- 90 40 48 152 37.6 9.3 Phe 150 330 281 39 259 141 254 Ant --- --- 34 25 135 33.7 9.8 Fla 500 740 583 101 599 353 8.3 Pyr 30 250 492 45 363 317 581 BaA 1620 2180 232 47 521 284 244 Cry 5400 3150 693 23 332 267 319 BbF 270 720 267 36 540 431 313 BkF 270 630 101 37 661 138 139 BaP 920 2030 136 35 461 294 249 Ind 250 750 47 --- 621 248 145 DaA 420 380 42 --- 623 120 21 Bghi 380 1150 70 --- 1618 221 239 Total 10270 12580 3077 885 7501 2834 3390 Table 2. Comparison of PAH concentrations in rural soils near to urban areas. Note: 1. This

Mexico, D.F.1 China2 India3, Delhi Korea4 Francia

zone. Highway Rural

zone

Park

Fig. 6. Distribution of PAHs in soils according to its number of aromatic rings in Tlahuac and Milpa Alta.

Important aspect for these soils in Tlahuac and Milpa are addition of farmyard manure and compost where we appreciated low values of PAHs in comparison crops (Cai et al., 2007). The organic matter plays a role to catch contaminants where can be stabilized within its structure. A long term, organic matter is long reservoirs of several contaminants which must be practice good management of this parameter of soil.

The presence of low molecular weight PAHs (2-3 aromatic rings) is indicative of petroleum (fossil fuel combustion) while high molecular weight PAHs (> 4 aromatic rings) are more likely to be derived from organic material combustion (Ma et al., 2005). There are likely to be multiple contamination sources in vegetables and garbage from these conservation terrains.

The contamination in Mexico City is high due to the environmental conditions, for example high population density and geography (such as altitude and temperature), where atmospheric deposition (diffuse and point sources) leads to high concentrations of PAHs. Similar contamination situations have been reported in Chinese soils (Ping et al., 2007).

In table 2, we appreciated high values of individual PAHs than others areas where compounds with high molecular weight surpass values of others studies of different countries. We supposed that high concentrations are due to pull of contaminants to South of City together certain activities of specific areas as possible explication.

Maliszewska-Kordybach (1996) suggested a soil contamination classification system based on ΣPAH16 as follows: noncontaminated soil (<200 ng/g), weakly contaminated soil (200– 600 ng/g), contaminated soil (600–1000 ng/g) and heavily contaminated soil (>1000 ng/g). According to this classification system, all samples in this study were heavily contaminated. This discrepancy among European and Mexican regulations are based in diverse criteria of analytical techniques mainly.

Studies made in past years for Mexico City has showed differences in abundance, type and distribution of PAHs according landscape type and territories of intensive urbanisation/

Milpa Alta Dry season

Fig. 6. Distribution of PAHs in soils according to its number of aromatic rings in Tlahuac

Important aspect for these soils in Tlahuac and Milpa are addition of farmyard manure and compost where we appreciated low values of PAHs in comparison crops (Cai et al., 2007). The organic matter plays a role to catch contaminants where can be stabilized within its structure. A long term, organic matter is long reservoirs of several contaminants which must

The presence of low molecular weight PAHs (2-3 aromatic rings) is indicative of petroleum (fossil fuel combustion) while high molecular weight PAHs (> 4 aromatic rings) are more likely to be derived from organic material combustion (Ma et al., 2005). There are likely to be multiple contamination sources in vegetables and garbage from these conservation terrains. The contamination in Mexico City is high due to the environmental conditions, for example high population density and geography (such as altitude and temperature), where atmospheric deposition (diffuse and point sources) leads to high concentrations of PAHs. Similar contamination situations have been reported in Chinese soils (Ping et al., 2007).

In table 2, we appreciated high values of individual PAHs than others areas where compounds with high molecular weight surpass values of others studies of different countries. We supposed that high concentrations are due to pull of contaminants to South of

Maliszewska-Kordybach (1996) suggested a soil contamination classification system based on ΣPAH16 as follows: noncontaminated soil (<200 ng/g), weakly contaminated soil (200– 600 ng/g), contaminated soil (600–1000 ng/g) and heavily contaminated soil (>1000 ng/g). According to this classification system, all samples in this study were heavily contaminated. This discrepancy among European and Mexican regulations are based in diverse criteria of

Studies made in past years for Mexico City has showed differences in abundance, type and distribution of PAHs according landscape type and territories of intensive urbanisation/

City together certain activities of specific areas as possible explication.

 Tlahuac Wet season

Milpa Alta Wet season

6 rings 5 rings 4 rings 3 rings 2 rings

0%

and Milpa Alta.

analytical techniques mainly.

 Tlahuac Dry season

be practice good management of this parameter of soil.

20%

40%

60%

80%

100%

industrialisation (Marr et al, 2004; Maliszewska-Kordybach et al., 2009). The spatial distribution of soils contaminated with PAHs reflected that rural/forested/ recreation areas decreasing along the South to North of Mexico City. Due in North there high concentration of industries and urban areas while in South of City the soils are considered conservation and recharge of aquifers with low density of population.

Table 2. Comparison of PAH concentrations in rural soils near to urban areas. Note: 1. This study, 2009; 2. Chen et al., 2005; 3. Agarwal, 2009; 4. Nam et al, 2003.

The importance of contamination of soil is due serious risks on population health due:

Further the contamination of water resource in form superficial or underground that are employ for irrigation for crops or human and animal consumption

Government hopes that in medium time the contaminants reduce significantly in favour the environment and health population of Mexico by programs of improvement of quality air.


**5.3.84 Trypan blue and trypan red**

The antimicrobial agents trypan blue and trypan red did not show a beneficial effect in ALS patients (Montanari & Pessina, 1955; Schwob & Bonduelle, 1952).

### **5.3.85 Vascular endothelial growth factor (VEGF, sNN0029)**

VEGF is a neuroprotectice and angiogenic growth factor (Maurer et al., 2008). It is currently tested in a phase I/II clinical trial (NCT00800501) for safety and tolerability. Of note, VEGF must be administered into the CSF.