LNG bunkering network design in inland waterways

Growing awareness of the environment and new regulations of the International Maritime Organization and the European Union are forcing ship-owners to reduce pollution. The use of liquefied natural gas (LNG) is one of the most promising options for achieving a reduction in pollution for inland shipping and short sea shipping. However, the infrastructure to facilitate the broad use of LNG is yet to be developed. We advance and analyze models that suggest LNG infrastructure development plans for refueling stations that support pipeline-to-ship and truck-toship bunkering, specifying locations, types, and capacities, and that take into account the characteristics of LNG, such as boil-off during storage and loading. We develop an effective primal heuristic, based on Lagrangian relaxation, for the solution of the models. We validate our approach by performing a computational study for the waterway network in the ArnhemNijmegen region in the West-European river network, including, among others, multi-year scenarios in which capacity expansion and reduction are possible.


Introduction
Environment-friendly transportation has become an important topic in the last decade. Many national and international policy makers encourage the use of alternative energy sources and fuels in the maritime sector. The European Commission promotes to develop the market for alternative transport fuels and the infrastructure to support them by the Alternative Fuels Directive 2014/94/ EU. This is in line with the Dutch energy agreement that sets an objective of reducing carbon dioxide emissions by 50% by 2030 (relative to markers in 2020) and in new international regulations for maritime transport (Burel et al., 2013;Arteconi and Polonara, 2013). Liquefied natural gas (LNG) is one of the most promising options for reducing pollution in the long term for short-sea and inland ships (SER, 2014). According to the European Union directive, by the end of 2025, countries must provide sufficient LNG stations at seaports, to accommodate LNG-powered ships.
The switch to the LNG fuel is also driven by the enforcement of the IMO (International Maritime Organization) emission standards. In October 2016, the leading regulatory body for the shipping industry's Marine Environment Protection Committee stated that the new sulfur emissions limits effective January 1, 2020 are 0.5 percent globally and 0.1 percent in IMO-designated emission control areas. The use of LNG instead of traditional Heavy Fuel Oil (HFO) reduces the emission of nitrogen oxides by approximately 80-85%, carbon dioxide by 20-30% (up to 40% for renewable LNG), and almost completely eliminates emissions of sulphur oxides. LNG as ship fuel is also able to compete with traditional fuel in areas beyond simple emissions levels: LNG has the potential to improve energy efficiency as well as to increase cost savings due to the lower operational costs of LNG shipsthe operating costs for a 33,000 DWT tanker ship can be reduced by 35% (Burel et al., 2013;Panteia, 2013). However, these cost savings only materialise after a certain payback period, depending on the cost of conversion from HFO to LNG (typically about 1.3 million Euros per ship), LNG price and infrastructure configuration (Burel et al., 2013). LNG-fuelled ships are more expensive to build, although costs continue to decrease as the technology becomes more mature. There are continued studies working towards developing systems for improved environmental footprint, a higher level of energy efficiency, as well as an improved boil-off rate and cargo capacity by considering actual operational conditions and optimisation in terms of hydrodynamics, machinery and system configuration (George et al., 2020;DNV GL, 2016a). The replacement of normal engines with the twice-as-expensive LNG engines in inland ships is already economically viable (Verbeek et al., 2011). Moreover, the price of LNG is expected to drop and become three times as cheap as crude oil within 30 years (DNV GL, 2016c). LNG-fuelled ships must adhere to certain safety regulations and standards along with the additional training and qualifications of the personnel involved (ISO, 2017). Accordingly, they are generally well designed, well maintained and operating with well-trained crew. Thus, LNG shipping so far has a good safety record (Vanem et al., 2008;Eliopoulou et al., 2016).
One drawback of LNG is the boil-off effect from storage and bunkering, which is the result of the low boiling point of natural gas. Another drawback is that, currently, an infrastructure of LNG bunker stations is not yet complete. This is the main reason that shipowners do not want to invest in LNG engines. On the other hand, entrepreneurs do not want to invest in LNG bunkering facilities, because at present there is little commercial demand for this type of fuel. This "chicken and egg" problem occurs with every new fuel introduction (Sperling and Kitamura, 1986). Consider, for example, the Dutch waterways. With 5000 km of waterway, almost 400 harbors, and one of the most modern inland fleets of 6900 vessels (Rijkswaterstaat, 2016), the market is ideally suited for LNG driven vessels. However, in 2017 only seven LNG fueled vessels were in use with plans for another 16 to be added in the near future (Port of Rotterdam, 2015;TNO, 2017), despite the goal of 100 LNG vessels and 500 LNG trucks by 2015 (DeltaLinqs, 2012). This slow progress is partly due to the lack of LNG-refueling stations along the Dutch waterways.
LNG can be provided to customers through different delivery options. Here, we consider the two main types of bunkering investments: pipeline-to-ship (PTS) and truck-to-ship (TTS). In pipeline-to-ship bunkering (also known as terminal-to-ship), LNG is bunkered directly from a terminal. Pipelines from the terminal to the vessel supply the LNG. In truck-to-ship bunkering, an LNG truck is connected to the receiving ship, using a flexible hose, assisted typically by a hose-handling manual cantilever crane (DNV GL, 2016b). TTS has lower infrastructure requirement and brings operational flexibility. However the capacity of TTS is limited and the supply of LNG for these trucks must be realized through a depot (Initiative WPC, 2015). The operating cost of truck refueling units is relatively low than the cost of operating a fixed facility.
In this environment, demand (for fuel) is characterized by a path from an origin to a destination, rather than a location. The goal is to deliver product along a predetermined route (Zeng et al., 2010). Sa and Ebrazi (2013) propose a formulation for a path (or flow) refueling location model, where the ranges of vehicles requiring the fuel are considered. In a setting closer to ours,  develop a model that maximizes the captured flow in a multi-period network by locating a predefined number of facilities for electric-car charging.
Previous models do not yet reflect many important characteristics of LNG waterway bunkering, e.g., boiling off of LNG during storage and loading and truck-to-ship bunkering. Truck-to-ship bunkering requires that a waterway network and a road network need to be considered simultaneously. The boil-off effect of LNG during storage and bunkering requires more careful modeling of inventory required to support fuel demand in the network. Furthermore, for strategic infrastructure planning, it is necessary to properly account for anticipated changes in fuel demand (which may be different for different corridors in the network considered).
In this research, we advance and analyze models that can assist government and companies in the design of an LNG distribution network. More specifically, we develop a multi-period capacitated demand capturing network design model considering pipeline to ship bunkering as well as truck-to-ship bunkering (TTS) and taking boil-off gas effects into account. To solve the models, we develop an effective primal heuristic, based on Lagrangian relaxation, which produces high-quality solutions, i.e., solutions that result in large profits (a high return on investment). We demonstrate the benefits of our model and solution approach by analyzing a portion of the West-European river network, the Nijmegen region, using real-life data, considering different scenarios and settings.
The remainder of the paper is organized as follows. Section 2 reviews the relevant literature. Section 3 introduces the mathematical model. Section 4 presents the Lagrangian based primal heuristic. Section 5 describes the results of a set of computational experiments to analyze the performance of the primal heuristic. Section 6 analyzes a section of the West-European River network. In Section 7, we present final remarks.

Literature review
Within academic research, a large number of articles have discussed facility location problems (Zvi Drezner, 2001). A comprehensive review of the existing literature in this area can be found in the article by Melo et al. (2009). Klose and Drexl (2005) categorized the literature into four broad categories: analytic models, continuous models, network models, and discrete models. A further distinction can be made based on how demand is structured, i.e., demand can be modelled by focusing on nodes or arcs. Traditional facility location models consider demand to be aggregated in nodes, and locate facilities on nodes. Another strategy is to relate the demand to arcs of the network, and assign the size of demand separately. These models identify paths between nodes of the network by balancing the demand coverage of the path with the cost of building sufficient fuel stations along the connecting arcs. The traditional facility location model uses population centers as demand sources, and enables travelling along an arc by locating the necessary number of fuel stations along that arc.
The demand structure analyzed in this research can be modelled with a path based setting. Demand is a traffic flow along a path covered by allocating a combination of facilities to the nodes along this path. These models are called flow capturing location models and were first used by Hodgson (1990). Flow is considered to be captured if a demand unit passes a facility (Schneider and Vis, 2016). The main reasoning behind the development of this type of model is that many types of facilities (e.g. ATMs and roadhouses) are not the main destination of trips (Berman et al., 1992;Zeng et al., 2009). These models are not yet suitable for refueling vehicles, because they assume that the flow is captured if facilities are situated anywhere along the path of a vehicle without considering the range of the vehicle. To address this, Kuby and Lim (2005), Yıldız et al. (2016), Upchurch and Kuby (2010) developed flow-refueling location models (FRLM) for road networks, which also considers the range of vehicles.
Research on FRLM is rather new and these studies focus on the location of alternative fuel (hydrogen, electric, and biofuel) facilities on a road network. The first FRLM study by Kuby and Lim (2005) required the pre-step of generating all valid refueling station combinations for each path, which results in a long solution time. Then, Capar and Kuby (2012) developed an efficient reformulation of the FRLM by introducing a decision variable that indicates if a vehicle at a certain node has enough remaining fuel to travel to the next open refueling station. This eliminates the generation of feasible facility combinations in the earlier model. Sa and Ebrazi (2013) proposed their own reformulation of the FRLM, which either determines a set of facility locations that minimizes the cost to serve a predetermined percentage demand (set covering) or maximizes the amount of captured flow or vehicle kilometre travelling for a given number of stations or a set budget (maximum coverage).
Next, Kim and Kuby (2012) proposed a variant of the FRLM where the model allows for deviation from the shortest path.  developed a multi-period model providing a forward-myopic and backward-myopic method as a solution technique. Recently, Kuby et al. (2017) used the standard FRLM model to support decision makers in building a fueling network for LNG trucks where country coverage as well as cross-border integration options are considered. Tran et al. (2018) developed a k-exchange neighbourhood heuristic to study the standard FRLM model.
The models above cannot be applied to our inland waterway setting. Accordingly we extend the current literature in the following aspects. We consider two main types of bunkering options: PTS and TTS. Our model considers a linked facility configuration in which trucks are replenished from the fixed bunkering facility. This features a new characteristic where inter-transshipment among different facilities are considered. We further consider capacity limitations in our model. The capacity of the linked facility must be sufficient for the supply of the truck facility. Here, two different transportation modes must be involved in the network as truck-to-ship bunkering uses the road network for transportation whereas the main demand sources from the waterways. The trucks are flexible and can be relocated based on shifting flow patterns. We also allow for the consideration of capacity decisions under the effects of boil-off separating the model used here from those that consider other fuels (e.g. hydrogen). The boil-off effect is influenced by the design of the tank and the composition of LNG (Querol et al., 2010).
The resulting model is challenging to be solved. We prove that our problem is NP-complete. We thereby develop a primal heuristic based on Lagrangian relaxation to efficiently and accurately solve the model.

Model description
We consider an inland waterway network D = (N, A) with node set N and arc set A. A node n ∈ N represents a possible location for a refueling station. An arc (i, j) ∈ A indicates the existence of a waterway (directly) connecting locations i and j. A set of shipping routes Q in the waterway network defines the demand. Each shipping route q ∈ Q, has an origin o q ∈ N, a destination d q ∈ N, a path p q = ((o q , and a demand f q , representing the number of ships expected to use the shipping route during the planning period. We assume that not all the demand has to be fulfilled. We further assume that a type z vessel consumes fuel at a rate v z (i, j) depending on the sailing speed on the arc (i,j), and that is refueled to its full capacity when it refuels. The latter simplifies the modeling, but also represents common practice. This allows us to define for each shipping route q ∈ Q, a network D q = (N q ,A q ), with N q the set of locations visited along path p q and A q a set of arcs (i, j) indicating that it is possible to reach j from i after refueling at i without refueling at intermediate locations. Fig. 1 shows an example of a shipping route from Düsseldorf to Worms in a small inland river network in which a vessel needs to be refuelled at Cologne and Mainz.
We consider two refueling options: (1) terminals providing PTS service, and (2) locations providing TTS service. The cost of providing PTS service at a location (i.e., operating and resupplying a terminal) is much higher than the cost of providing TTS service at a location (i.e., transporting product by truck from a terminal to that location). This makes TTS service an interesting option as it can be used to expand coverage to locations with smaller demands at a relatively low price. We let L⊆N be the set of locations where it is possible to build a terminal and provide PTS service.
We consider evaporation, as it is a critical aspect in the supply of LNG, and, therefore, has to be modeled in order to provide a more accurate representation of LNG inventories and LNG deliveries. More specifically, we capture the boil-off effect during storage, loading, and transportation with rates ψ, ϕ and η, respectively. For ease of presentation, we introduce the model assuming a single period and a full tank at the origin. Thereafter, we extend to a multi-period model.
To model the possibility of only partially serving the demand of shipping route q, we modify D q and introduce a dummy source s and sink k and arcs (s, o q ), (d q , k) and (s,k). Next, we let the variable x q ij (0⩽x q ij ⩽1) represent the fraction of the demand of shipping route q that uses arc (i, j) ∈ A q . If x q sk > 0, we interpret the value of x q sk as the fraction of the demand that is not served. Furthermore, we let binary variable y ic indicate whether a terminal with capacity c ∈ C is operated at location i ∈ L, where C is the set of possible terminal capacities.
To model the possibility that trucks transfer LNG from one location to another, we introduce the variables w ij for i ∈ L and j ∈ N, indicating the mass transferred from a terminal operated at location i to a TTS station at location j during the planning period. The supply capacity of the terminal at location i should be sufficient to satisfy the requirement of the TTS station at location j.
We assume a limited budget is available to cover the operating cost during the planning period. We let I denote the total budget and we let m denote the maximum fraction of the budget that can be allocated to TTS operations.
The parameters and decision variables are summarized below (relative to the planning period): Sets: Set of nodes associated with shipping route q A q Set of arcs associated with shipping route q C Set of facility types (i.e., set of capacities) L Set of nodes where it is possible to operate a terminal Z Set of vessel types Parameters: Distance between locations i and j in the water network (km) fq Demand on shipping route q (number of vessels) v z (i, j) LNG consumption (tons/km) of vessel of type z on the trip from node i to j p Selling price of LNG (per ton) ψ Boiling off rate while in storage (% per ton) ϕ Boiling off rate while loading (% per ton) η Boiling off rate while in transport (% per ton) a(c) Product supply capacity of a facility of type c (tons) u(i, c) Cost of operating and resupplying a facility of type c at location i Cost of transporting product from location i to location j (per ton) I Budget m Maximum fraction of the budget that can be allocated to TTS operations

Decision variables:
yic Binary variable indicating whether a terminal with capacity c is operated at location The fraction of the demand of shipping route q that uses arc (i, j) wij Amount of LNG (in tons) transported from a terminal at location i to a TTS station at location j The model can now be stated as: The objective function maximizes profit and has three terms. The first term captures the revenue from LNG sales. Because we assume that a ship's tank is always filled up to capacity, it suffices to know the previous location where a ship refueled to calculate the amount of LNG consumed, and, thus, the revenue; the quantity loaded is equal to the distance traveled multiplied by the fuel consumption rate (adjusted by the fraction of the demand served). The second term captures the cost of operating and resupplying the terminals that we decide to operate. The third term captures the cost of supplying locations providing TTS service using over-the-road truck transport. Constraints (2) represent flow conservation. Constraints (3) ensure that no more LNG is loaded at a location than is available. The left-hand side of the inequality represents the total amount of LNG loaded at the location taking into account losses due to boil-off, i.e., it captures the fact that to load a quantity W, a quantity 1 1− ϕ W of the available supply is consumed (due to boiling-off during loading). Note that the boiling off rate used in our model is an average number. The right-hand side of the inequality represents the availability of LNG at the location taking into account (a) boiling-off during storage (in case of PTS service) and truck transportation (in case of TTS service) and (b) product taken from (in case of PTS service) and delivered to (in case of TTS service) the location. Constraints (4) ensure that locations providing TTS service only receive product from locations providing PTS service. Constraint (5) and (6) capture the budget constraints.
The model presented above is intended to be used as part of a larger investment decision process and focuses on analyzing potential profits as a function of the operating budget. As such, it ignores all aspects related to capital expenditures and, for example, does not restrict the number of terminals to be opened. Enforcing a condition on the number of terminals to be opened, e.g., a specific number or a maximum number, can easily be incorporated if desirable.
The following theorem shows that the problem is NP-complete.

Theorem 1. The single-path LNG network design problem is NP-complete.
Proof. From PARTITION, i.e., given a set of positive integers b 1 ,b 2 ,…,b n , with ∑ n i=1 b i /2 = η, we check whether there exists a subset S⊆ {1, …, n} such that ∑ i∈S b i = η. Let the single path along which vessels travel have n potential locations for refueling stations, and let the journey from the start of the path to the end of the path be such that it requires (at least) one refueling stop for a vessel, but that fueling stop can be at any of the n potential locations. Furthermore, let the capacity of a facility operational at location i ∈ {1, …, n} be such that b i vessels can be served, let the cost of operating a facility at location i ∈ {1, …, n} be ρb i with 0 < ρ < 1, and let the budget be ρη. Finally, let the demand, i.e., the number of vessels traversing the path, be η and the revenue from the sale of fuel to a vessel traversing the path be 1. Then the maximum profit is achieved when all vessels are served, which means that fueling capacity for η vessels has to be available along the path, which is possible if and only if PARTITION has a feasible solution. □ Further, we discuss two possible extensions of the above model. First, the single period model can be extended to a multi-period model with T planning periods by adding a time index to the variables, and adding the following set of additional constraints: Constraints (9) enforce that terminals cannot suspend operations once operations have started. Second, we extend our model to allow for modular capacities. We allow to adjust capacities and therefore eliminate constraint (9) and redefine y ic as y c ′ ic and u(i, c) as u c ′ (i,c). y c ′ ic equals 1 if the capacity of the facility located at i is changed from type c to type c ′ where c, c ′ ∈ C. We define C = {0,C}where 0 means there is no operating facility. u c ′ (i, c) is the cost of changing the capacity at location i from type c to type c ′ .

A primal heuristic
Solving the formulation presented in Section 3 with realistic size instances can be (too) time-consuming. Therefore, in this section, we focus on developing a primal heuristic, based on Lagrangian relaxation, which can produce high-quality solutions in a reasonable amount of time. We want to clarify that the choice of the initial feasible solution of the primal heuristic is designed based on the unique characteristics of the model followed by a Lagrangian procedure where we update the primal bound by solutions generated by the relaxation.
We first provide a high-level overview of the solution approach and then continue with explaining the key elements in more detail. A high-level overview of the solution approach is presented in Fig. 2 and contains the following steps: 1. Construct an initial primal solution (Section 4.1). 2. Initialize the Lagrangian multipliers and set the stopping criteria (Section 4.2). 3. Solve the Lagrangian relaxation to obtain a dual bound. 4. Obtain a primal solution by solving the original formulation with an additional constraint that restricts the choice of facilities to open to those that had a positive value in the solution to the Langrangean relaxation. Update the primal bound if an improved solution is found. 5. Terminate if the one or more of the stopping criteria are met, and continue if not. 6. Update the Lagrangian multipliers using either the Polyak or the heavy ball step size rule and go to Step 3.
Next, we discuss the node-to-path heuristic and the key elements of the solution approach in more detail.

An initial feasible solution
An initial feasible solution is created by greedily selecting locations to open facilities so as to capture as-yet uncovered demand as long as any budget remains. For each selected location, we assign a facility with the highest capacity. We consider two schemes for calculating a location's potential.
The first scheme is based on the following observations: (1) a long shipping route brings in more revenue, (2) opening a facility at a location is more attractive if a large number of shipping routes pass (and may visit) that location, and (3) opening a facility at a location is more attractive if is more likely that the shipping routes that pass the location can be "captured", i.e., if no, or only a few, additional facilities need to be opened along the shipping route. For location i ∈ L, let Q i ⊆Q denote the set of shipping routes that pass location i.
For each shipping route q ∈ Q i , let N q denote the set of locations along the route that do not yet have a facility. The potential of a This gives the following heuristic: 1. P ← ∅; 2. For i ∈ L⧹P, determine ϒ i , the location potential; 3. i * ← argmax i∈U⧹P ϒ i ; 4. P ← P ∪ {i * }; Fig. 2. High-level overview of the primal heuristic.

E. Ursavas et al.
Transportation Research Part C 120 (2020) 102779 7 5. If all shipping routes are captured, then stop. If not and there is remaining budget, then go to Step 2.
The second scheme is based on the following observation. A shipping route is attractive if the revenue it can generate is high and the cost to capture the shipping route is low. The revenue a shipping route can generate is determined by the number of ships, the length of the route, and the price of fuel. The cost to capture a shipping route is determined by the additional facilities that need to be opened along the route to make it possible for the ships to go from origin to destination. The latter can be found by solving a shortest-path problem (in which location with already opened facilities can be visited for free). This gives the following heuristic: 2. For q ∈ Q, determine S q the minimum cost set of locations that need to be opened to capture q (taking into account the locations and associated facilities in P that have already been selected); 5. If all shipping routes are captured, then stop. If not and there is remaining budget, then go to Step 2.

A Lagrangian relaxation
Lagrangian relaxation is a decomposition technique in which a set of (complicating) constraints is relaxed, but violations of these constraints are penalized in the objective function. The coefficients in the objective function that represent these penalties are also known as the Lagrange multipliers. The Lagrangian dual problem is to find the set of multiplier values that gives the tightest bound. See Fisher (2004) for an introduction to Lagrangian relaxation.
When applying Lagrangian relaxation it is critical to select an appropriate set of constraints to relax (to dualize). The selection has to balance the quality of the resulting bound and the computation time needed to solve the relaxed problem. A preliminary analysis and investigation has shown that dualizing the constraints (9) provides a good compromise between the quality of the bound and the computational effort required to compute the bound. The objective value of this relaxed problem will be a dual bound on the optimal value of the original problem. In computing the Lagrangian multipliers, we use the Polyak and heavy-ball techniques (Polyak, 1987).
In a subgradient method with a Polyak update, Lagrange multipliers γ (the gradient) for constraints (9) are updated as follows: where θ is calculated as with Z is the objective value of the best feasible solution found so far, Z iter the objective value of the relaxed problem, and δ (scalar between 0 and 2) the agility parameter. With a heavy-ball update, Lagrange multipliers γ (the gradient) for constraints (9) are updated as follows: where α and s are calculated as with 0⩽β⩽1 the filter parameter (which controls the amount of memory used by the algorithm). All Lagrange multipliers are initialized to zero.

Numerical experiments
We first demonstrate the performance of our primal heuristic and then conduct a scenario analysis of the West-European river network. All computational experiments are performed using a 64-bit Windows server 2008 R2 operating system on an Intel Xeon CPU E5-2660 2.2 GHz processor with 8 Gb of RAM. Computation times reported in the results tables are in seconds. In the result tables "Compact" indicates that the model presented in Section 3 was solved using a commercial solver whereas "LB-B" and "LB-P" indicate that the primal heuristic was used with heavy-ball and Polyak updates, respectively. An initial primal solution is constructed using the two heuristics given in Section 4 (selecting the best among the two solutions). The maximum number of iterations in the primal heuristic is set to 30. For all the experiments in this section, we consider two types of vessels. The speed of a vessel on the arc (i, j) is generated randomly based a uniform distribution in the range from 10 to 30.
Figs. 3 and 4 illustrate the behavior of the heuristic for an instance with 30 nodes, 100 arcs, and 6 time periods for both Polyak and heavy-ball updates. We show the value of the lower and upper bound in each iteration (as well as the best known lower and best known upper bound at each iteration). We see that for update schemes, the heuristic converges after 15 iterations. The gap between the dual bound (UB) and the primal bound (LB) is less than 0.1 percent. The decrease in the dual bound is significant already after two iterations. The behavior for other instances show a similar convergence pattern. This demonstrates that the heuristic is effective.
Next, we analyze the effect of the network size on the performance of the different solution approaches. The number of nodes and arcs is varied for a setting with 15 different origin-destination pairs and 8 time periods. The number of nodes per path in each of these instances varies between five and nine and the number of paths visiting a node varies between two and five. The results can be found in Table 1, where the column labeled #N/P reports the average number of nodes in demand paths and the column labeled #P/N reports the average number of demand paths visiting a node. We observe that the time required to solve the mixed integer programming formulation increases rapidly as the size of the network increases and the number of nodes in a demand path increases. We also see an increase in the solution time of the heuristic, but it is much smaller. The quality of the solutions produced by the heuristic deteriorates, but it is still capable of producing high-quality solutions.
Next, considering the initial scenario with 25 nodes, 70 arcs, 15 origin-destination pairs and 8 time periods as a base case, we test the computational performance for different number of arcs. In practice a network with a higher number of arcs corresponds to a geographical setting with an equal number of ports and possible facility locations, but a higher number of rivers connecting to them. These results showing the computational performance with the increase in the number of arcs can be found in Table 2.
We continue by exploring the effect of the number of time periods on the performance of the different solution approaches. The results can be found in Table 3. We see that a small increase in the number of time periods leads to a notable increase in solve time for the mixed integer programming formulation, whereas it has little effect on the solution time of the heuristics and on the quality of the solutions they produce.
Next, we test the model with the increasing number of O-D pairs. In Table 4 the computational results are shown for a network with 25 nodes, 70 arcs, and 8 time periods with different number of O-D pairs. We observe the increase in the computational time required for increasing number of O-D pairs for the compact model, similar to the previous experiments. Both primal heuristics are robust to these changes.
All in all, we observe that the number of time periods, network size, and number of O-D pairs affect the solving time of the compact model significantly. Compared to the compact model the primal heuristics perform better in computation time. For large networks and/or numerous time periods, both primal heuristics are recommended due to the good performance in terms of both computation speed and accuracy. Both the techniques find solutions with a gap of at most 5 %.

Case study: West-European River network
We consider the West-European river network and conduct scenario analysis to evaluate current and future LNG bunkering supply network design choices. The Nether Rhine and Waal are two of the most important rivers in the European waterway network. A significant amount of traffic between the sea ports and the German inland ports along the Rhine as well as from Basel in Switzerland uses these rivers. Fig. 5 shows a map of the West-European river network including cities where PTS and TTS services exist and/or are The characteristics of the area, such as the distances through the waterway and road networks and the demand in the area, have been obtained from the Expertise-en Innovatie Centrum Binnenvaart (EICB), the Municipality of Nijmegen, and through interviews with shipping companies operating in these waterways as part of the thesis study by Gollackner (2015). In order to obtain the required inputs for our model, a combination of empirical and qualitative data collection methods have been used. To determine the potential demand for LNG, inland waterway traffic has been analyzed to identify vessels suitable for conversion to LNG. The data collected included ship Maritime Mobile Service Identity, latitude and longitudes, speed, course, status. The traffic data passing the Arnhem-    Nijmegen region was then matched with detailed information about the individual ships, e.g., the bunkering volume of the ship, the year the ship was built, and the ship's dimensions, in order to be able to assess the feasibility of converting to LNG. Furthermore, interviews with stakeholders have been conducted to gather information about behavior and preferences of ship owners, and to identify region-specific constraints for our distribution network model. Four shipping companies have been studied in detail to further gain an understanding of the potential demand and requirements for an LNG distribution network. These shipping companies operate barges that pass through the Arnhem-Nijmegen region. We only consider barges with a minimum length of 135 meters and a remaining lifetime of more than 20 years, which would be suitable for conversion to LNG. The traffic patterns of these barges were used to derive the set of paths to be covered. This case study has 163 shipping services (O-D pairs), 43 locations, and 44 links between locations. We consider a planning horizon with 10 periods, each period representing a year. The demand is assumed to be stable across the periods. Small budget increases occur in each period. We consider three types of PTS stations, each with different capacities. The cost of operating a PTS station depends on its location and the facility type. The cost of operating a TTS station depends of the PTS station it is connected to (capturing the cost of driving trucks back and forth between the PTS station and the TTS station). Due to the size of the instance, it was only possible to obtain solutions using the heuristic presented in Section 4. For each computational experiment, we present a table specifying the locations selected and the facilities operated at these locations. For example, Rheinberg(3) indicates that a PTS station with a facility of Type 3 is operated in location Rheinberg, and 5 Duisburg [Rheinberg] indicates that a TTS station is operated in location Duisburg where 5 units of LNG is supplied from Rheinberg. There are already operational LNG bunkering stations in the Port of Rotterdam, Amsterdam, and Antwerp, so we have assumed these in our model, i.e., our focus is on how to extend the existing LNG supply network.
We start by presenting the LNG supply network suggested by our model for West-European River network. We follow that by an analysis of the benefits of incorporating TTS service in the supply network. Next, we study the effect of available budget, and, finally, we analyze a situation in which capacity adjustments are possible (i.e., we assume modular PTS facilities) to accommodate expected changes in the demand pattern due to the dredging of the River Scheldt. Table 5 shows the proposed supply network for the West-European River for periods 3, 6 and 10, in which major changes to the supply network occur. Total profit is 24 million euros. We observe TTS stations are intensively used: 28 TTS stations are deployed in period 10. Furthermore, we see that two primary regions are served. Initially (until period 6), the area between Rheinberg and Mainz is served, but afterwards also the Dutch part of the network, the area between Amsterdam and Vlissingen is served. Fig. 6 provides a graphical representation of the LNG supply network in the final period (arc labels indicate the sailing distance between nodes).
Next we study the benefits of incorporating TTS stations. To do so, we eliminate the budget constraint and compare the solution in which the possibility of TTS stations is included to the one in which the possibility of TTS stations is not included. Table 6 summarizes the characteristics of the networks in the final period. The top part of the table shows the characteristics of the network when TTS stations are an option and the bottom part of the table shows the characteristics of the network when TTS stations are not an option. When the budget is not limiting, we have a 41.7 percent increase in the profits in the network with TTS stations capturing approximately 70 percent of the total demand. In this setting, the ports of Rheinberg, Mainz, Dusseldorf, Cologne, and Vlissingen are chosen with large (Type 3) facilities (in addition to the existing facilities at the Port of Rotterdam, Amsterdam, and Antwerp). From these facilities, 38 TTS stations are supplied. This suggests that operating TTS stations is cheaper than operating a PTS station with a small (Type 1) facility. An example is the TTS station in Leverkussen, which is not far from Cologne, where 7 units of LNG is supplied from Cologne. Likewise, in Rotterdam we observe that the present PTS station is not sufficient and we further need to increase the supply. This is realized with TTS supplied from Vlissingen. Without the option to open TTS stations, PTS facilities are chosen in Rheinberg (3), Rotterdam(2), Krefeld(1), Cologne(2), and Mainz (3). Not only are the chosen locations different, but smaller facilities are operated in these locations. Importantly, fewer shipping demands are captured (10 percent fewer), resulting in lower profit (approximately 8 percent less). For example, the shipping demand between Antwerp and Duisburg is not served. The results also show the need for further extending the current facilities at the Port of Rotterdam. Clearly, the benefit of incorporating TTS facilities in the LNG network for the West-European River network are significant.
Next, we investigate the impact of the budget (i.e., the yearly budget increase) on profit. More specifically, we vary the yearly budget increase and observe the changes in profit. The results can be found in Fig. 7. We observe that the increasing rate in profit drops by about a half when the yearly budget increase exceeds 700,000 Euros.
Finally, we investigate the West-European River network when accounting for the planned dredging operations in River Scheldt. With the increased depth resulting from dredging, it is expected that the number of vessels visiting the Antwerp harbor will increase as the ships will have a smoother passage with fewer delays due to depth and tide conditions. At the same time, competing ports may experience a decrease in demand. These changes will require capacity adjustments at the bunkering stations. To allow for these adjustments, we slightly modify our model and analyze the LNG supply network when the anticipated changes in demand are incorporated. More specifically, after the dredging of the Scheldt, the Port of Antwerp will be able to serve larger vessels and some demand will likely shift towards the Port of Antwerp from the Port of Rotterdam (currently the only harbor that can accommodate large vessels). Our results show that the increase in demand in the Antwerp region has a considerable impact on the LNG supply network. We see that LNG facilities are added in Dordrecht and Moerdijk to capture shipping routes via Antwerp. Furthermore, we see a decrease in the capacity required in Rotterdam.

Final remarks
Use of LNG, as a marine fuel is expected to increase. Therefore, a significant growth in the development of LNG supply network is anticipated. In line with this expectation, this paper studied the design of LNG network in an inland waterway network motivated from initiations towards a more environmentally friendly water transportation in the West-European rivers.
Existing studies on network design were not applicable for the LNG maritime environment since they did not reflect many important features of LNG waterway bunkering. Therefore we proposed a mixed integer programming model that features LNG characteristics such as terminal-to-ship, truck-to-ship bunkering and boil-off gas effects. An effective primal heuristic based on Lagrangian relaxation is developed that provides solutions that result in large profits. We tested the performance of our solution technique on different river networks and demand patterns such as networks with a geographically dense canals, higher number of ports and different shipping services.
In Section 6 we conducted a case study based on the inland waterway traffic in the West-European river network using data of the stakeholders of the LNG network. The model forms the required LNG supply network in this region with regard to the demand along the waterway paths. Next to the harbors at the Port of Rotterdam, Amsterdam and Antwerp, additional LNG stations are suggested mainly in Rheinberg, Cologne, Mainz, Duisburg and Vlissingen with differing capacities under restricted budget. These locations along with the truck facilities connected to them are promising in terms of providing higher profits. We further analyzed the possible outcomes after the dredging of River Scheldt. In such a scenario we expect increase in required facility capacities in Dordrecht and Moerdijk whereas a reduction at the facility capacity in Rotterdam would be profitable to eliminate higher operating costs. Our focus has been on analyzing the potential profits of an LNG bunkering network in inland waterways given an operating budget for a planning period. This involves locating and sizing bunkering stations as well as deciding on the type of bunkering station, i.e.,  Fig. 6. Facility placement plan for the West-European river network.

Table 6
Budget unrestricted outcome PTS TTS pipeline-to-ship and truck-to-ship. The results of a case study in the West-European region show that the use of truck-to-ship stations in addition to pipeline-to-ship stations is highly beneficial. We conclude with a few observations. Clearly, understanding the potential profit of an LNG bunkering network given an operating budget for a planning period is only part of the investment decision process. The construction of LNG terminals requires large investments to cover the costs of acquiring permits, site preparation, foundations, piping, terminal building, etc. Furthermore, in practice, LNG services are likely to represent only a portion of a terminal's operations. Investments in LNG bunkering networks may not happen if government do not come to the table with tax exemptions and subsidies. This makes it necessary to consider the direct and indirect socio-economic impact LNG facilities bring to an area. This shows that the development of new LNG infrastructures should be assessed comprehensively which involves more than just economic considerations. According to Council of European Energy Regulators (CEER) to foster growth of the LNG market, increasing LNG flexibility and liquidity, removing rigid contracts clauses, simplifying Third Party Access (TPA) to gas storage conditions, suggesting new services to foster terminal utilisation and exploring ways to promote price signals are utmost important.
(CEER, 2019). The potential profit of an LNG bunkering network given an operating budget for a planning period provides a critical piece of information in the overall investment decision process. Total profit 10 11 Fig. 7. Impact of total budget on total profit.