Carbon Dioxide Storage Site Location and Transport Assignment Optimization for Sustainable Maritime Transport

Abstract

Maritime carbon dioxide (CO₂) transport plays a pivotal role in facilitating carbon capture and storage (CCS) systems by connecting emission sources with appropriate storage sites. This process often incurs significant transportation costs, which must be carefully balanced against penalties for untransported CO₂ resulting from cost-driven decisions. This study addresses the CO₂ storage site location and transport assignment (CSSL-TA) problem, aiming to minimize total tactical costs, including storage site construction, ship chartering, transportation, and penalties for direct CO₂ emissions. We formulate the problem as a mixed-integer programming (MIP) model and demonstrate that the objective function exhibits submodularity, reflecting diminishing returns in facility investment and ship operations. A case study demonstrates the model’s effectiveness and practical value, revealing that optimal storage siting, strategic ship chartering, route allocation, and efficient transportation significantly reduce both transportation costs and emissions. To enhance practical applicability, a two-stage planning framework is proposed, where the first stage selects storage sites, and the second employs a genetic algorithm (GA) for transport assignment. The GA-based solution achieves a total cost only 2.4% higher than the exact MIP model while reducing computational time by 57.9%. This study provides a practical framework for maritime CO₂ transport planning, contributing to cost-effective and sustainable CCS deployment.

1. Introduction

With the accelerating global climate crisis, the need to achieve deep decarbonization across energy and industrial systems has become a central concern of governments and industries. As widely acknowledged by the Intergovernmental Panel on Climate Change and the International Energy Agency, carbon capture and storage (CCS) is indispensable for achieving net-zero emission targets, especially in hard-to-abate sectors such as cement, steel, and chemicals [1,2,3]. In this context, the construction of reliable carbon dioxide (CO₂) storage and transport infrastructure has emerged as a critical enabler of large-scale CCS deployment. Particularly in Europe, where emission sources are scattered across the continent while geological storage capacity is concentrated in offshore basins such as the North Sea, efficient planning of CO₂ shipping infrastructure is vital [4].

Maritime CO₂ transport is increasingly recognized as a vital component in carbon capture and storage (CCS) due to its inherent flexibility and cost efficiency [5]. Compared to fixed pipelines, ship-based transport offers significant advantages, particularly during the initial phases of CCS deployment and in regions with spatial and temporal variability [6]. Maritime transport generally requires lower initial capital investment than pipelines, which demand substantial upfront costs, making ships more suitable for projects with uncertain demand evolution [7]. Pipelines are generally more efficient for shorter distances with high CO₂ volumes [8], with Knoope et al. [7] specifically showing that pipelines are preferred for distances under 200 km when dealing with CO₂ volumes greater than 5 million tons/year. Conversely, for smaller volumes, such as 2.5 million tons/year, shipping becomes competitive at distances greater than 500 km. Furthermore, ships offer operational flexibility, often requiring less frequent maintenance than pipelines, which need regular inspections, especially in corrosive environments [7]. Unlike pipelines, which may require costly structural reinforcements to withstand harsh environmental conditions, ships provide adjustable routing to adapt to operational needs [6]. This combination of reduced capital commitment, operational scalability, adaptability to variable conditions, and suitability for long-distance transport makes maritime transport a compelling and flexible alternative to pipelines for linking industrial clusters to offshore storage hubs [5]. However, the economic viability of ship transport largely depends on the strategic deployment of coastal storage stations, which serve as essential docking and injection points in the maritime supply chain. Inappropriate placement of these facilities may result in extended detours, under-utilized fleet capacity, increased berthing fees, and even emission penalties due to CO₂ that remains untransported [9,10]. Thus, the integrated planning of storage station locations, fleet deployment, and shipping schedules is of central importance for CCS operational efficiency and policy compliance.

Beyond economic and strategic considerations, maritime CO₂ transport also faces distinct thermophysical and operational constraints that significantly affect system design and deployment. CCS is pivotal for mitigating climate change by capturing CO₂ from industrial sources and storing it in geological formations like deep saline aquifers or depleted hydrocarbon fields, supporting global net-zero emission goals [11]. Maritime transport is essential for connecting emission sources to offshore storage sites, especially for long distances (>500 km) or lower CO₂ volumes (<2.5 million tons/year), where it outperforms pipelines in cost effectiveness [12]. However, CO₂ maritime transport requires maintaining specific pressure (0.7–1.5 MPa) and temperature (223–246 K) conditions to keep CO₂ in a liquid state, ensuring high density (≈1159 kg/m³) and avoiding solidification or supercritical transitions that could cause operational issues like blockages [13,14]. These physical and engineering constraints significantly influence storage site selection, as sites often require supercritical CO₂ for injection, and ship-route allocation, which must account for vessel design and dynamic factors like weather or port congestion [15]. The intricate interplay of these thermophysical requirements and operational considerations presents substantial challenges for the efficient and safe design of CCS networks. The necessity to meticulously manage CO₂ physical state across the supply chain, from transport to injection, while considering logistical complexities, underscores the critical need for robust optimization approaches. Therefore, developing frameworks for the strategic layout of storage sites and transport routes that explicitly account for these factors is paramount to enhance the safety, efficiency, and cost effectiveness of CCS supply chains, thereby facilitating their global deployment. This study is motivated by these challenges, proposing an optimization framework that specifically integrates CO₂ physical state requirements into strategic decision-making for site selection and route allocation.

To complement these technical and operational challenges, several real-world CCUS projects provide practical validation and experience that reinforce the necessity of integrated site selection and transport planning. Saudi Aramco’s Uthmaniyah CO₂ Europe Demonstration Project and Norway’s Northern Lights Project are two prominent examples. Since 2015, the Uthmaniyah project has captured approximately 800,000 tonnes of CO₂ annually from the Hawiyah Natural Gas Liquids plant, compressing it to 1500–1600 psi and transporting it via an 85 km pipeline to a Jurassic carbonate reservoir for injection [16]. Advanced monitoring ensures secure storage, with about 40% of CO₂ permanently sequestered [17]. The project’s success reflects effective infrastructure design, cross-disciplinary integration, and commitment to environmental goals [18], offering transferable insights into offshore site planning and pipeline-based CO₂ transport. In contrast, the Northern Lights Project, operational since 2024, adopts a fully maritime approach: liquefied CO₂ is shipped 100 km from Øygarden to a North Sea saline aquifer using 7500 m³ vessels, with a total capacity of 1.5 million tonnes per year [19]. Monitoring via seismic and pressure sensing ensures long-term storage security [20]. The project’s ship-based flexibility, cross-border coordination, and policy support from Norway’s Longship initiative [21] provide a working reference for marine-based CCS logistics. Systematically identifying the factors contributing to the success of these projects reveals several dimensions that are directly relevant to the proposed optimization model. These include the importance of strategic site selection based on storage capacity and geographic efficiency, the allocation of flexible ship-based transport routes to accommodate industrial demand, and the role of operational constraints such as berthing frequency and infrastructure reliability. The Northern Lights project demonstrates the feasibility of ship-route planning under variable conditions, which aligns with the model’s capability to assign ships, optimize transport frequency, and support multi-source routing. Uthmaniyah’s focus on reservoir suitability and monitoring emphasizes the value of matching infrastructure investment with storage performance, a core aspect of the site selection sub-model. Together, these projects inform both the structural formulation and the practical applicability of the proposed approach. Their design and implementation experiences directly parallel the optimization problems explored in this study, particularly regarding routing strategies and offshore injection planning.

Nevertheless, determining optimal locations for CO₂ storage sites is challenging due to the spatial heterogeneity of emission sources, port infrastructure, and long-distance maritime routes. Existing studies focus extensively on the optimization of CO₂ transportation supply chains using pipelines and ships [22,23,24]. These studies explore various aspects of supply chain design, including the integration of capture, transport, and storage processes, the selection of efficient transport modes, and the development of resilient networks to support CCS deployment. Additionally, research has investigated large-scale infrastructure network design at national or regional levels [25], compared the economic and logistical feasibility of pipelines versus ships [26], and conducted techno-economic analyses to assess the cost effectiveness and environmental impact of ship-based CO₂ transport [27,28]. Some studies have also developed generalized routing strategies for broad geographical areas [29]. However, few studies integrate the strategic siting of coastal CO₂ storage stations with the optimization of detailed maritime transport operations, such as ship scheduling, task assignment, and route planning, to efficiently connect multiple emission sources to these sites. For instance, some frameworks optimize CCS networks by selecting storage sites but prioritize pipeline-based systems, with limited attention to maritime logistics like ship scheduling or multi-source collection [30]. Other studies analyze ship-based CO₂ transport, addressing port infrastructure and ship design, but often neglect critical operational constraints such as the rational allocation of transportation routes based on the volume of CO₂ generated [27,31].

Although these studies offer important insights, they typically overlook key operational constraints that are critical in maritime CO₂ transport planning. First, they often fail to model the required berthing frequency for each emission source, which is essential to ensure that annually generated CO₂ can be transported without delays or accumulation. Second, most models assume full transport or no transport, whereas in real applications, a portion of CO₂ may remain untransported and incur emission penalties, which must be incorporated explicitly in cost evaluation. Third, the existing routing models rarely allow a single ship to collect CO₂ from multiple sources (i.e., dual-source routing), despite its operational relevance in reducing underutilization and improving routing flexibility. These operational features can lead to unrealistic or infeasible plans, highlighting a research gap that this study aims to address.

To address these gaps comprehensively, this study formulates a novel CO₂ storage site location and transport assignment (CSSL-TA) problem. The CSSL-TA problem jointly optimizes the strategic construction of CO₂ storage sites within coastal maritime transport systems and the operational assignment of ship routes and transport tasks across multiple emission sources. A mixed-integer programming (MIP) model is developed to support integrated decision-making, simultaneously determining which storage sites to construct, how to assign ships to feasible transport routes, and how to allocate annual CO₂ volumes from each emission source under system-wide constraints. These realistic constraints include ship transport capacity limits, storage site capacity constraints, berthing frequency requirements at each emission port, annual route operating limits, and CO₂ direct emission constraints. These constraints ensure that the solution is both operationally feasible and aligned with emission coverage regulations. Overall, the proposed MIP model provides a practical and implementable optimization framework for supporting long-term infrastructure planning in coastal CO₂ transport networks.

In particular, the multifaceted contributions of this paper are outlined as follows:

We establish an MIP model for jointly optimizing the location of CO₂ storage sites, the assignment of ships to feasible transport routes, and the annual number of round trips made by each chartered ship. The realistic operational constraints in real-world CCS operations are well incorporated in our model. This formulation enables stakeholders to minimize the total cost while supporting practical infrastructure planning for long-term maritime CO₂ transport systems.

We prove the necessity of maintaining integrality for construction, ship assignment, and routing frequency decisions. We also demonstrate that the objective function is submodular with respect to all three variable types. This finding highlights the inherent diminishing returns effect in the CO₂ transport system, whereby additional station construction, ship-route assignments, or round trips contribute progressively less to overall cost reduction.

Sensitivity analyses on annual available sailing time, sailing speed, and the capacities of candidate storage sites are conducted to examine their effects on total cost, facility deployment, and key performance indicators. The results reveal several important insights: a longer annual sailing time allows ships to complete more round trips and reduces the number of ships required; faster sailing speeds increase transport efficiency but may raise fuel costs and direct emissions; and larger storage capacity reduces the number of required storage sites but may lead to longer transport distances due to wider spatial distribution. These findings provide practical guidance for balancing transport efficiency and infrastructure layout under different operational assumptions.

A case study demonstrates the practical value of the proposed model, offering actionable managerial insights into storage construction, ship-route assignment, and operational planning. This study shows that (i) selecting appropriate storage site locations helps reduce transport distances and lower transportation costs; (ii) strategically chartering different types of ships and assigning them to suitable routes effectively reduces direct CO₂ emissions; (iii) ensuring that ships operate the appropriate number of times on selected routes guarantees the timely transportation of CO₂ from each emission source, and helps minimize total cost by balancing CO₂ transport cost with direct emission penalties.

To improve the model’s applicability in industries, we further propose a two-stage planning framework that separates strategic storage site selection from transport assignment. The first stage determines site construction decisions by balancing the supply–demand relationship between storage capacities and CO₂ generation volumes, alongside construction costs at candidate sites, while the second stage applies a genetic algorithm (GA) to solve the transport assignment problem. This framework offers a practical alternative to the exact MIP model and can serve as a useful decision-support tool in the scenarios where computational efficiency is critical.

The remainder of this paper is structured as follows. Section 2 reviews the related work on CCS supply chain optimization and CO₂ transport. Section 3 formulates the research problem as an MIP model and presents its theoretical properties. Section 4 reports numerical results and sensitivity experiments. Section 5 introduces the two-stage optimization framework and details the GA-based solution approach. Finally, Section 6 concludes this paper and outlines directions for future research.

2. Literature Review

To comprehensively understand the methodologies for CO₂ storage site selection and transport assignment in CCS systems, this review examines the literature focusing on optimization models and algorithms for CSSL-TA within maritime transportation frameworks.

2.1. CO₂ Storage Site Selection

The selection of geological storage sites for CO₂, particularly deep saline aquifers and depleted hydrocarbon fields, typically utilizes multi-criteria decision-making frameworks combined with thorough geological assessments. Bachu [31] establishes fundamental geological screening criteria, including sufficient reservoir depths to ensure CO₂ remains in a dense state, adequate porosity and permeability, and the presence of a competent caprock with sufficient thickness and low permeability to ensure containment. Grataloup et al. [32] develop a two-stage screening approach employing both exclusionary and qualification criteria, utilizing Geographic-Information-System-based spatial analysis for site identification in the Paris Basin. Alcalde et al. [33] propose a four-stage workflow integrating Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) analysis to systematically weigh criteria such as storage capacity, injectivity potential, and containment risk within the Acorn project context. Raza et al. [34] emphasize the importance of petrophysical models, like Archie’s equation for water saturation determination and Bahadori’s correlation for predicting CO₂ density under reservoir conditions, alongside understanding various trapping mechanisms relevant to saline aquifers. Hsu et al. [35] apply the Analytic Network Process to prioritize interdependent criteria, notably caprock permeability, thereby enhancing the robustness of site evaluation. For offshore contexts, Luo et al. [36] demonstrate the use of seismic inversion and well-log data integration to assess structural trap integrity and reservoir-caprock characteristics in the South Yellow Sea Basin, combining geophysical and petrophysical analyses. Following International Organization for Standardization guidelines, Aarnes et al. [37] employ probabilistic risk assessments and dynamic modeling to evaluate potential CO₂ migration pathways and seal integrity. Furthermore, Akhurst et al. [38] describe risk-assessment-led pre-characterization methodologies that utilize publicly available data in the early stages to identify and mitigate potential risks associated with storage security, thereby guiding subsequent technical characterization efforts. Techno-economic assessments by Gruson et al. [39] on four diverse European storage sites confirm that storage costs are highly site-specific, influenced significantly by the chosen development strategy and monitoring plan, highlighting the difficulty in establishing a generic storage cost. Table 1 summarizes key studies on CO₂ storage site selection, highlighting their considered factors and applicable scopes.

2.2. Optimization of Maritime CO₂ Transport

The optimization of maritime CO₂ transport uses advanced modeling techniques to navigate logistical complexities. Techno-economic modeling by Roussanaly et al. [12], optimizing ship parameters and port infrastructure over significant distances, indicates that maritime transport presents advantages over pipelines for lower CO₂ volumes. Zhang et al. [9] develop an MIP model to optimize the strategic deployment of a CCS supply chain in Northeast China. This model integrates the selection of emission sources, capture technologies, pipeline routes, intermediate hubs, utilization options such as Enhanced Oil Recovery, and geological storage sites, demonstrating potential net cost reduction through revenue generation. Addressing the need for robust planning under uncertainty, Elahi et al. [40] formulate a multi-stage stochastic MIP optimization framework for the UK CCS system. This approach accommodates temporal and spatial integration while considering financial market and storage capacity uncertainties, yielding flexible investment strategies designed to minimize potential financial losses. Knoope et al. [41] employ Real Options Analysis (ROA) to evaluate the impact of uncertainties like CO₂ pricing and transportation tariffs on infrastructure investment timing, finding that such uncertainties can significantly elevate the required pre-investment price threshold. To enhance supply chain resilience, Gabrielli et al. [42] propose an MIP model that minimizes expected unsequestered CO₂ under various disruption scenarios. Utilizing scenario-based stochastic programming, their model incorporates alternative transport options and demonstrates potential cost savings relative to deterministic planning. Engebø et al. [43] introduce a risk-based optimization framework integrating fault tree analysis and Monte Carlo simulations to quantitatively assess disruption risks and inform the prioritization of redundant shipping routes. Focusing on large-scale network development, Kjärstad et al. [44] integrate techno-economic modeling of the European power sector with detailed CO₂ transport infrastructure analysis. Their work investigates potential CCS network configurations spanning the whole of Europe, emphasizes how the spatio-temporal distribution of capture facilities and site-specific storage properties crucially determine optimal pipeline routing and deployment schedules, and quantifies the significant cost impact of restricting onshore storage options. Similarly, Morbee et al. [45] detail the InfraCCS model, a tool designed using methodological innovations such as node clustering and specialized algorithms for route pre-selection, to determine the optimal least-cost extent and associated investment cost of a potential large-scale European CO₂ pipeline transport network required for decarbonization scenarios up to 2050. Employing methodological innovations for route selection and costing, the model estimates network investment requirements, primarily for pipelines, and highlights the critical need for cross-border coordination. Targeting the specifics of maritime logistics, Bennæs et al. [46] formulate an MIP-based model to optimize key parameters of a ship-based CO₂ transportation system, including fleet composition, interim storage, and ship specifications. Implemented using Gurobi, this model minimizes overall logistics costs while adhering to operational constraints like port berthing availability and delivery schedules. Complementing transport optimization, Ajayi et al. [47] review the methodologies for estimating CO₂ storage capacity, and advocate for the integration of both static and dynamic models within multi-criteria decision analysis frameworks. Incorporating Monte Carlo simulations to manage geological uncertainties, their approach demonstrates potential improvements in site ranking accuracy. Table 2 summarizes key studies on optimization of maritime CO₂ transport, highlighting their objectives, methodologies, and limitations.

Our work distinguishes itself from existing research by adopting a more comprehensive approach to optimizing the CO₂ chain. As highlighted in the comparative table, our model is unique in its simultaneous consideration of four key objective costs: CO₂ transport cost, CO₂ direct emission cost, ship-chartering cost, and storage site construction cost. In contrast, previous studies typically focus on a subset of these factors. For instance, while many studies address CO₂ transport and storage site construction costs, none of the cited works incorporate all four critical cost components that our research addresses. This holistic consideration, facilitated by our modeling method employing MIP with a two-stage optimization using a genetic algorithm, allows for a more robust and realistic assessment of the overall CCTS project economics and environmental impact.

2.3. Research Gap

Despite these significant advances, several key gaps remain in the existing literature. First, the majority of the studies do not adequately address critical operational constraints specific to maritime transport, such as limitations on berthing frequency, the explicit handling of penalties for direct CO₂ emissions, and the complexities of multi-source, multi-destination routing. Second, there is a lack of integrated models that simultaneously optimize strategic site selection and detailed transport assignment within a single framework. Third, efficient solution methods are often not proposed, especially for tackling the computational challenges posed by large-scale instances. To bridge these gaps, we propose an MIP model designed for the CSSL-TA problem, which integrates strategic site selection with detailed transport assignment. This MIP model uniquely incorporates real-world operational constraints related to berthing frequency, direct emission penalties, and multi-source routing to achieve a minimal-cost and operationally feasible solution. Additionally, recognizing large-scale computational demands, we introduce a practical two-stage solution framework where a GA efficiently optimizes transport assignment following initial site selection, thereby facilitating more effective planning for sustainable CCS deployment. These all confirm the statement of our contributions to the existing literature described at the end of Section 1.

3. Problem Formulation

In this section, we first introduce the problem background and describe the challenges in CO₂ storage site location and ship-routing assignment in Section 3.1. Following that, we formulate the problem using an MIP model in Section 3.2. Finally, Section 3.3 provides a detailed description of model analysis.

3.1. Problem Description

The escalating need to mitigate CO₂ emissions has positioned CCS as a cornerstone technology in the global transition to a low-carbon economy. Governments and industries alike are increasingly driven by stringent environmental regulations and ambitious climate targets, spurring significant investments in CCS systems. These systems not only capture CO₂ from industrial processes but also prevent its release into the atmosphere, thereby addressing both environmental and regulatory challenges. The development of CCS infrastructure entails a coordinated interplay of technological innovation, financial investment, and policy support. Such integration is critical to the creation of a robust, scalable CCS network that can operate efficiently across long-term planning horizons and manage large quantities of CO₂ in a cost-effective and environmentally sustainable manner.

In maritime CCS systems, the strategic placement of coastal CO₂ storage sites presents a complex optimization challenge driven by several competing factors. Constructing these sites requires significant investment. Meanwhile, ships transporting CO₂ operate along predefined routes that link emission sources to selected storage sites. If storage sites are poorly positioned, ships may need to follow longer routes, leading to increased fuel consumption. These inefficiencies elevate transportation costs and may also result in higher direct CO₂ emissions, which incur substantial penalty costs under emission control regulations. In light of these considerations, it is crucial to precisely determine the optimal locations of CO₂ storage sites and ship assignments to minimize the total system cost in the planning horizon. This cost includes the construction cost of selected CO₂ storage sites, the chartering cost of ships (incurred when a ship is chartered for CO₂ transportation service), the transportation cost (including fuel cost along sailing routes and berthing cost of ships), and the penalty cost for direct CO₂ emissions into the atmosphere. Meanwhile, the system must satisfy annual CO₂ capture and emission constraints: the total CO₂ emissions generated at each emission source in a year must be transported to a storage site in a timely manner; otherwise, some of them will be released to the atmosphere within the permitted limit.

Consider the CO₂ storage site location problem in a coastal-area transportation network during a multi-year planning horizon, where the number of years in the planning horizon is denoted by $T$. The set of candidate locations where CO₂ storage sites can be constructed is denoted by $I$ (indexed by $i$). Hereafter, we use ‘storage site’ and ‘location’ interchangeably. The set of CO₂ emission sources is denoted by $P$ (indexed by $p$). Ships are used to transport CO₂ from the emission sources to storage sites along fixed transport routes. The set of all feasible CO₂ transport routes is denoted by $R$ (indexed by $r$). Each route $r\in R$ includes exactly one storage site $i\in I$, and either one or two emission sources from $P$, i.e., $i\left(r\right)\in I$, $P\left(r\right)\subseteq P$, with $\left|P\left(r\right)\right|\in \left\{1, 2\right\}$. If a route connects one emission source with a storage site, it represents a round trip between the two; if it connects two emission sources with a storage site, it represents a round trip in which the ship departs from the storage site, visits both emission sources in sequence, and returns to the storage site.

The set of ships available for charter is denoted by $S$ (indexed by $s$). The decision of whether to charter a ship is based on the actual transport demand. If a ship is chartered, it remains in service for the entire planning horizon, and its chartering cost $C_s^{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{p}}$ is incurred as a one-time fixed expense. Once a ship is chartered, it is assigned to a single route $r\in R$, and this assignment remains fixed in the planning horizon. Each ship is required to berth both when receiving CO₂ at an emission source and when injecting CO₂ to a storage site. The berthing cost of ship $s$ on route $r$, denoted by $B_{sr}$, depends on the number of emission sources visited along the route. Each berthing operation of ship $s$ at an emission source or a storage site incurs a fixed cost $b_s$, which includes tugboat assistance, pilotage, and mooring services. If route $r$ includes two emission sources, then $B_{sr}=3b_s$; if it includes only one emission source, then $B_{sr}={2b}_s$. This parameter contributes to the transportation cost together with the fuel cost $F_{sr}$. For a route involving one emission source and one storage site, a full round trip consists of traveling from the emission source to the storage site and back. For a route that includes two emission sources and one storage site, a full round trip consists of sailing from the storage site to the first emission source, continuing to the second emission source, and then returning to the storage site. The maximum number of round trips that ship $s$ can make on route $r$ per year is denoted by $N_{sr}$, which depends on the ship’s sailing speed $v_s$ and the route’s sailing distance $l_r$. Assuming a full year is available for operation and denoting one year’s available sailing time by $t_0$, $N_{sr}$ is calculated as the available sailing time in one year divided by the time required for one round trip on route $r$, i.e., $N_{sr}=t_0/\left(l_r/v_s\right)$. This parameter is not required to be an integer, as any incomplete trip at the end of a year is assumed to be completed in the following year. To account for realistic factors that significantly impact maritime CO₂ transport, such as weather conditions, route congestion, and berthing delays, it is important to recognize their influence on ship sailing speeds and available navigation time. Adverse weather, congestion, and delays can reduce sailing speeds or decrease the available navigation time and limiting the $N_{sr}$. Adjusting $N_{sr}$ to reflect these operational uncertainties and constraints can effectively capture the real-world challenges faced in maritime shipping.

Each storage site $i\in I$ has a capacity $V_i$, which represents the maximum amount of CO₂ that can be injected into the site per year. The construction cost of a candidate storage site $i$, denoted by $C_i$, consists of two components: a fixed cost component, incurred only if the location $i$ is eligible for construction, and a variable cost component $b_i·V_i$, which scales linearly with the site’s storage capacity, i.e., $C_i=a_i+b_i{V}_i$. The amount of CO₂ produced by each emission source $p\in P$ is fixed per year and denoted by $E_p$. To ensure timely removal of CO₂ from each source, a minimum berthing frequency requirement $N_p$ is imposed, specifying the minimum number of ship arrivals required at source $p$ annually. For each round trip, the amount of CO₂ carried from source $p$ by ship $s$ must not exceed the ship’s transport capacity $Q_s$. If any CO₂ remains untransported by the end of the year, it is directly released into the atmosphere. The amount of CO₂ that emission source $p$ is allowed to release per year is capped by $E_p^{\mathrm{m}\mathrm{a}\mathrm{x}}$, and the penalty cost of emitting a ton of CO₂ into the atmosphere is denoted by $\widehat{C}$.

The objective function of this study is to minimize the total system cost in the planning horizon, including the construction cost of the selected CO₂ storage sites, the total chartering cost of all chartered ships, the total transportation costs consisting of both fuel and berthing expenses, and the penalty costs for CO₂ that is directly emitted into the atmosphere. To this end, we need to decide which candidate storage sites to construct. This decision is represented by the binary variable $x_i$, which equals 1 if the $i$-th candidate CO₂ storage site is constructed, and 0 otherwise. We also determine whether each ship $s\in S$ is chartered and, if so, which route it is assigned to. This is represented by the binary decision variable $z_{sr}$, which equals 1 if ship $s$ is chartered and assigned to route $r$, and 0 otherwise. The number of round trips that ship $s$ makes on route $r$ per year is described by the continuous decision variable $n_{sr}$. The amount of CO₂ transported from source $p$ on route $r$ by ship $s$ per year is denoted by the continuous decision variable $w_{srp}$. The amount of CO₂ released into the atmosphere by source $p$ per year is represented by the continuous decision variable $e_p$.

3.2. Model Design

In this subsection, we present an MIP model based on the problem setting described above. The model captures the integrated planning decisions involved in the construction of CO₂ storage sites and the chartering and routing of ships for CO₂ transport. It focuses on strategic decisions for constructing storage sites and operational decisions concerning ship assignment, route scheduling, and annual CO₂ transportation volumes. The objective is to minimize the total system costs while adhering to practical constraints such as storage capacity limits, transport capabilities requirements, berthing frequency requirements, and emission regulations. Table 3 summarizes the notations used in the model.

The MIP model is developed to optimize the construction of CO₂ storage sites, the assignment of ships to transport routes, and the annual routing frequency of each ship, which can be written as follows:

min	${\displaystyle \sum _{i\in I}}{C}_i{x}_i+{\displaystyle \sum _{s\in \mathrm{S}}}{\displaystyle \sum _{r\in R}}{C}_s^{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{p}}{z}_{sr}+{\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r\in R}}T\left(F_{sr}+B_{sr}\right)n_{sr}{z}_{sr}+{\displaystyle \sum _{p\in P}}T\widehat{C}{e}_p$		(1)
s.t.	${\displaystyle \sum _{r\in R}}{z}_{sr}\le 1$	$\forall s\in S$	(2)
	${\displaystyle \sum _{p\in P\left(r\right)}}{w}_{srp}\le Q_s{n}_{sr}{z}_{sr}$	$\forall s\in S, \forall r\in R$	(3)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:i\left(r\right)=i}}{\displaystyle \sum _{p\in P\left(r\right)}}{w}_{srp}\le V_i{x}_i$	$\forall i\in I$	(4)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:p\in P\left(r\right)}}{n}_{sr}{z}_{sr}\ge N_p$	$\forall p\in P$	(5)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:p\in P\left(r\right)}}{w}_{srp}+e_p=E_p$	$\forall p\in P$	(6)
	$0\le n_{sr}\le N_{sr}{z}_{sr}$	$\forall s\in S, \forall r\in R$	(7)
	$x_i\in \left\{0,1\right\}$	$\forall i\in I$	(8)
	$z_{sr}\in \left\{0,1\right\}$	$\forall s\in S,\forall r\in R$	(9)
	$w_{srp}\ge 0$	$\forall s\in S, \forall r\in R, \forall p\in P$	(10)
	$0\le e_p\le E_p^{\mathrm{m}\mathrm{a}\mathrm{x}}$	$\forall p\in P$	(11)

The objective function (1) minimizes the total cost, which consists of four components: the construction cost of selected CO₂ storage sites, the total chartering cost of all chartered ships, the transportation cost in the planning horizon (including fuel and berthing costs), and the penalty cost for direct CO₂ emissions into the atmosphere. Constraint (2) ensures that each ship can be assigned to at most one route, reflecting the condition that not all ships must be chartered. Constraint (3) denotes that for each ship assigned to a route, the total amount of CO₂ transported per year does not exceed the product of its transport capacity and the number of annual round trips. Constraint (4) requires that the total injected CO₂ into each constructed storage site per year does not exceed its capacity. Constraint (5) requires that each emission source is visited by ships a sufficient number of times per year to meet its minimum berthing requirement, thereby ensuring an adequate transport frequency and enabling timely removal of the generated CO₂. Constraint (6) ensures that all CO₂ generated by each source is either transported or emitted. Constraint (7) imposes operational bounds on the annual number of round trips that each chartered ship can make on its assigned route, ensuring that this value does not exceed the maximum allowable trips based on ship speed and route distance. Constraints (8)–(11) are the domains of decision variables.

3.3. Model Analysis

In this section, we introduce some theoretical properties regarding the MIP model presented in Section 3.2.

3.3.1. Linearization of Bilinear Terms Involving Ship Assignment Variables

To simplify the structure of the MIP model and reduce computational complexity, we next examine whether the bilinear term $n_{sr}{z}_{sr}$ can be safely linearized. Through analytical derivation, we establish Lemma 1.

Lemma 1.

The bilinear terms$n_{sr}{z}_{sr}$in the model can be replaced by linear terms$n_{sr}$.

Proof. 

In the MIP model, the bilinear term $n_{sr}{z}_{sr}$ appears in both the objective function to calculate total cost and Constraints (3) and (5) to limit the total transported volume by the ship’s capacity and to ensure the minimum berthing requirements for each source by controlling the number of annual round trips. According to Constraint (7), $0\le n_{sr}\le N_{sr}{z}_{sr}$, the number of round trips $n_{sr}$ is bounded by the binary variable $z_{sr}$. This directly implies that when $z_{sr}=0$, $n_{sr}=0$; and when $z_{sr}=1$, $n_{sr}\in \left[0, N_{sr}\right]$. In both cases, the product $n_{sr}{z}_{sr}$ is is exactly equal to $n_{sr}$. Therefore, according to the feasible region defined by Constraint (7), we have: $n_{sr}{z}_{sr}=n_{sr}$, $\forall s\in S, \forall r\in R$, which completes the proof. □

This replacement eliminates unnecessary nonlinearity from both the objective and the constraints, improving computational tractability while preserving all model logic.

3.3.2. Submodular Property of the Objective Function

A set function $f:2^Y\to \mathbb{R}$ defined on the subsets of a finite set $Y$ is called submodular if it satisfies the diminishing returns property: for every $Y_1\subseteq Y_2\subseteq Y$ and every $y\in Y\setminus Y_2$, it holds that $f(Y_1\cup \{y\left\}\right)-f\left(Y_1\right)\ge f(Y_2\cup \{y\left\}\right)-f\left(Y_2\right)$ [48]. In simple terms, a submodular function is a set function that describes the relationship between a set of inputs and an output, where adding an additional input has a decreasing additional benefit (diminishing returns). This property essentially means that the marginal gain from adding an element to a smaller set is at least as great as the gain from adding the same element to a larger set. Submodular functions are particularly interesting because they help model a variety of naturally occurring phenomena, such as economies of scale, network effects, or the spread of information in social networks.

For notation convenience, we denote $\mathit{x}={\left(x_i\right)}_{i\in I}$, $\mathit{z}={\left(z_{sr}\right)}_{ s\in S, r\in R}$, $\mathit{n}={\left(n_{sr}\right)}_{ s\in S, r\in R}$ and $\mathit{e}={\left(e_p\right)}_{ p\in P}$, which are restricted in the feasible region $\mathsf{\Omega }$. Then, the objective function can be presented by a function of $\mathit{x}$, $\mathit{z}$, $\mathit{n}$ and $\mathit{e}$, denoted by $F_{\Omega }(\mathit{x},\mathit{z},\mathit{n},\mathit{e})$. We then give the following lemma about $F_{\Omega }(\mathit{x},\mathit{z},\mathit{n},\mathit{e})$. The feasible region Ω is defined by Constraints (2)–(11), which include operational and capacity constraints such as construction decisions, ship chartering and route assignments, ship round-trip limits, emission allocations, and berthing frequency requirements. The objective function $F_{\Omega }(\mathit{x},\mathit{z},\mathit{n},\mathit{e})$ represents the total cost in the planning horizon, composed of storage site construction cost, ship chartering cost, transportation cost (including fuel and berthing cost), and penalty cost for untransported CO₂ emissions, where the decision variables are constrained in the feasible region Ω.

Lemma 2.

The objective function${ F}_{\Omega }(\mathit{x},\mathit{z},\mathit{n},\mathit{e})$is submodular regarding variables$\mathit{x}, \mathit{z}$and$\mathit{n}$.

Proof. 

The objective is to minimize the total cost of selecting a subset of locations $x_i$ for storage site construction, a subset of ships and route assignments $z_{sr}$, and a subset of round-trip quantities $n_{sr}$, to cover all required service elements and satisfy the direct emission limits. The direct emission volume $e_p$ is linked to the transport CO₂ volume $w_{srp}$ via Constraint (6): ${\sum }_{s\in S}{\sum }_{r:p\in P\left(r\right)}{w}_{srp}+e_p=E_p$, $\forall p\in P$. The set covering function is known to be submodular [48], which suggests that our objective function $F_{\Omega }(\mathit{x},\mathit{z},\mathit{n},\mathit{e})$ is also submodular with respect to the variables x, z, and n.

In the subsequent contents, we prove this submodularity by a more rigorous step-by-step analysis.

Step 1: The binary variables $x_i$ indicate whether the CO₂ storage site at location $i$ is constructed. Define ${\mathit{x}}_1$ and ${\mathit{x}}_2$ as two subsets of constructed CO₂ storage sites and let ${\mathit{x}}_1\subseteq {\mathit{x}}_2\subseteq \mathit{x}$, where each $x_i\in \mathit{x},\forall i\in I,$ is 1. Then, for any $x_{i^{\prime }}\notin {\mathit{x}}_2$, the marginal gain by adding another storage site $x_{i^{\prime }}$ to ${\mathit{x}}_1$ can be defined as $F_{\mathsf{\Omega }}\left({\mathit{x}}_1,\mathit{z},\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left({\mathit{x}}_1\bigcup \left\{1\right\},\mathit{z},\mathit{n},\mathit{e}\right)$. We note that the marginal gain is represented by $F_{\mathsf{\Omega }}\left({\mathit{x}}_1,\mathit{z},\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left({\mathit{x}}_1\bigcup \left\{1\right\},\mathit{z},\mathit{n},\mathit{e}\right)$ since the objective function is the cost. Similarly, $F_{\mathsf{\Omega }}\left({\mathit{x}}_2,\mathit{z},\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left({\mathit{x}}_2\bigcup \left\{1\right\},\mathit{z},\mathit{n},\mathit{e}\right)$ is the marginal gain by adding a same storage site.

Step 2: Since ${\mathit{x}}_1\subseteq {\mathit{x}}_2$, any potential cost reduction in service coverage by adding $x_{i^{\prime }}$ to ${\mathit{x}}_2$ is also possible when added to ${\mathit{x}}_1$, but may remain unchanged due to already covered demands with the same cost. Therefore, the marginal gain by adding another storage site $x_{i^{\prime }}$ to the existing subsets of sites ${\mathit{x}}_1$ is greater than or equal to adding it to ${\mathit{x}}_2$.

Step 3: In the presence of limited injection demands, the possibility of yielding zero marginal benefit by adding a new storage site $x_{i^{\prime }}$ to ${\mathit{x}}_2$ is higher than adding it to ${\mathit{x}}_1$. In other words, when most emission sources are already well served in ${\mathit{x}}_2$, adding a new site will provide negligible additional benefit, thus reinforcing the diminishing return property.

Conclusively, we always have $F_{\mathsf{\Omega }}\left({\mathit{x}}_1\mathit{z},\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left({\mathit{x}}_1\bigcup \left\{1\right\},\mathit{z},\mathit{n},\mathit{e}\right)\ge F_{\mathsf{\Omega }}\left({\mathit{x}}_2,\mathit{z},\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left({\mathit{x}}_2\bigcup \left\{1\right\},\mathit{z},\mathit{n},\mathit{e}\right)$.

Similar submodular behavior holds for the variables $z_{sr}$ and $n_{sr}$, as shown in the following analysis.

The binary variables $z_{sr}$ indicate whether ship $s$ is chartered and assigned to route $r$. Define ${\mathit{z}}_1$ and ${\mathit{z}}_2$ as two subsets of ship-route assignments and let ${\mathit{z}}_1\subseteq {\mathit{z}}_2\subseteq \mathit{z}$, where each $z_{sr}\in \mathit{z},\forall s\in S, \forall r\in R$ is 1. Then, for any $z_{{sr}^{\prime }}\notin {\mathit{z}}_2$, the marginal gain by adding another ship-route assignment $z_{{sr}^{\prime }}$ to ${\mathit{z}}_1$ can be defined as $F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_1,\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_1\bigcup \left\{1\right\},\mathit{n},\mathit{e}\right)$. We note that the marginal gain is represented by $F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_1,\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_1\bigcup \left\{1\right\},\mathit{n},\mathit{e}\right)$ since the objective function is the cost. Similarly, $F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_2,\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_2\bigcup \left\{1\right\},\mathit{n},\mathit{e}\right)$ is the marginal gain by adding the same storage site. Since ${\mathit{z}}_1\subseteq {\mathit{z}}_2$, any potential cost reduction in service coverage by adding $z_{{sr}^{\prime }}$ to ${\mathit{z}}_2$ is also possible when added to ${\mathit{z}}_1$, but may remain unchanged due to already covered demands with the same cost. Therefore, the marginal gain by adding another ship-route assignment $z_{{sr}^{\prime }}$ to the existing subsets of sites ${\mathit{z}}_1$ is greater than or equal to adding it to ${\mathit{z}}_2$. In the presence of limited transport demands, the possibility of yielding zero marginal benefit by adding a new assignment $z_{{sr}^{\prime }}$ to ${\mathit{z}}_2$ is higher than adding it to ${\mathit{z}}_1$. In other words, when most emission sources are already well served in ${\mathit{z}}_2$, adding a new site will provide negligible additional benefit, thus reinforcing the diminishing return property. Conclusively, we always have $F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_1,\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_1\bigcup \left\{1\right\},\mathit{n},\mathit{e}\right)\ge F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_2,\mathit{n},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},{\mathit{z}}_2\bigcup \left\{1\right\},\mathit{n},\mathit{e}\right)$.

Additionally, the continuous variables $n_{sr}$ indicate the number of round trips that ship $s$ make on route $r$. Define ${\mathit{n}}^{\left(1\right)}$ and ${\mathit{n}}^{\left(2\right)}$ as two feasible round-trip decisions and let ${\mathit{n}}^{\left(1\right)}\le {\mathit{n}}^{\left(2\right)}$. Then, the marginal gain by increasing a round trip quantity $n_{{sr}^{\prime }}$ from the lower level ${\mathit{n}}^{\left(1\right)}$ can be defined as $F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(1\right)},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(1\right)}+n_{{sr}^{\prime }},\mathit{e}\right)$. Similarly, $F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(2\right)},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(2\right)}+n_{{sr}^{\prime }},\mathit{e}\right)$ is the marginal gain by increasing a round trip from the higher level ${\mathit{n}}^{\left(2\right)}$. Since ${\mathit{n}}^{\left(1\right)}\le {\mathit{n}}^{\left(2\right)}$, any potential cost reduction brought by increasing the number of round trips from a lower level ${\mathit{n}}^{\left(1\right)}$ is also possible when increasing from the higher level ${\mathit{n}}^{\left(2\right)}$ but may remain unchanged due to already fulfilled transportation tasks. Therefore, the marginal gain by increasing the number of round trips from ${\mathit{n}}^{\left(1\right)}$ to ${\mathit{n}}^{\left(1\right)}+n_{{sr}^{\prime }}$ is greater than or equal to increasing it from ${\mathit{n}}^{\left(2\right)}$. In the presence of limited transportable volumes, the possibility of yielding zero marginal benefit by increasing a new round trip quantity $n_{{sr}^{\prime }}$ from ${\mathit{n}}^{\left(1\right)}$ is higher than increasing it from ${\mathit{n}}^{\left(2\right)}$. In other words, when most transport volumes are already well served by a large number of round trips ${\mathit{n}}^{\left(2\right)}$, adding more offers limited additional benefit, thus reinforcing the diminishing return property. Conclusively, we always have $F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(1\right)},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(1\right)}+n_{{sr}^{\prime }},\mathit{e}\right)\ge F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(2\right)},\mathit{e}\right)-F_{\mathsf{\Omega }}\left(\mathit{x},\mathit{z},{\mathit{n}}^{\left(2\right)}+n_{{sr}^{\prime }},\mathit{e}\right)$.

The above analysis verifies the diminishing return property for the decision variables $x_i$, $z_{sr}$ and $n_{sr}$. Therefore, the objective function is submodular regarding these variables. □

3.3.3. Model Tightening via Capacity–Demand Coupling

To further strengthen the LP relaxation of the model and reduce the solution space without excluding any feasible integer solution, we examine whether constraints involving capacity limits can be tightened based on the actual amount of CO₂ allocated in the system. By exploiting the property of demand–capacity coupling, we establish Lemma 3 and Lemma 4.

Lemma 3.

The right-hand side of the annual storage capacity Constraint (4) can be tightened by replacing${ V}_i$with$\min \left(V_i,E_i\right)$, where$E_i$represents the total amount of CO₂ generated by all emission sources connected to the storage site$i$.

Proof. 

Consider Constraint (4) of the model: ${\sum }_{s\in S}{\sum }_{r:i\left(r\right)=i}{\sum }_{p\in P\left(r\right)}{w}_{srp}\le V_i{x}_i$, $\forall i\in I$, where $w_{srp}$ indicating the amount of CO₂ transported from source $p$ on route $r$ by ship $s$ per year, and $V_i$ is the capacity of the $i$ -th CO₂ storage site. We denote the set of emission sources that are connected to storage site $i$ by $P\left(i\right)$. Specifically, $P\left(i\right)\subseteq P$. The total amount of CO₂ generated by all emission sources in $P\left(i\right)$ is expressed as follows: $E_i={\sum }_{p\in P\left(i\right)}{E}_p$. From Constraint (6), we know that ${\sum }_{s\in S}{\sum }_{r:p\in P\left(r\right)}{w}_{srp}+e_p=E_p$, $\forall p\in P$, which ensure that all CO₂ generated by each source is either transported or emitted. Therefore, for each constructed site $i$ ($x_i=1$), the value $E_i$ equals the sum of the CO₂ volume transported from emission sources $p\in P\left(i\right)$ to all storage sites connected to these sources, plus the total CO₂ volume directly emitted from $p\in P\left(i\right)$, that is:

(12)$\sum _{s\in S}\sum _{r\in R}\sum _{p\in P\left(i\right)}{w}_{srp}+\sum _{p\in P\left(i\right)}{e}_p=E_i.$

Since each emission source connects to at least one storage site, the following inequality holds:

(13)$\sum _{s\in S}\sum _{r:i\left(r\right)=i}\sum _{p\in P\left(r\right)}{w}_{srp}\le \sum _{s\in S}\sum _{r\in R}\sum _{p\in P\left(i\right)}{w}_{srp}.$

Combined with (12), we have:

(14)$\sum _{s\in S}\sum _{r:i\left(r\right)=i}\sum _{p\in P\left(r\right)}{w}_{srp}+\sum _{p\in P\left(i\right)}{e}_p\le E_i.$

According to Constraint (11), $e_p$ is non-negative for all $p\in P$, and thus we can conclude that ${\sum }_{r:i\left(r\right)=i}{\sum }_{p\in P\left(r\right)}{e}_p\ge 0$. Then, we have ${\sum }_{s\in S}{\sum }_{r:i\left(r\right)=i}{\sum }_{p\in P\left(r\right)}{w}_{srp}\le E_i{x}_i$ for $x_i=1$ and ${\sum }_{s\in S}{\sum }_{r:i\left(r\right)=i}{\sum }_{p\in P\left(r\right)}{w}_{srp}=0$ for $x_i=0$. Conclusively, we always have ${\sum }_{s\in S}{\sum }_{r:i\left(r\right)=i}{\sum }_{p\in P\left(r\right)}{w}_{srp}\le E_i{x}_i$, $\forall i\in I$. Therefore, we can replace $V_i$ in the original Constraint (4) with $\mathrm{m}\mathrm{i}\mathrm{n}\left(V_i,E_i\right)$. The modified constraints become

(15)$\sum _{s\in S}\sum _{r:i\left(r\right)=i}\sum _{p\in P\left(r\right)}{w}_{srp}\le \mathrm{m}\mathrm{i}\mathrm{n}\left(V_i,E_i\right)·x_i, \forall i\in I,$

which completes the proof. □

Lemma 4.

The right-hand side of the ship transport capacity Constraint (3) can be tightened by replacing$Q_s{n}_{sr}$with$\mathrm{m}\mathrm{i}\mathrm{n}\left(Q_s{n}_{sr}, E_r\right)$, where$E_r$represents the amount of CO₂ generated by all emission sources on route$$.

Proof. 

Consider Constraint (3) of the model: ${\sum }_{p\in P\left(r\right)}{w}_{srp}\le Q_s{n}_{sr}{z}_{sr}$, $\forall s\in S, \forall r\in R$, where $w_{srp}$ indicates the amount of CO₂ transported from source $p$ on route $r$ by ship $s$ per year, and $Q_s$ is the transport capacity of ship $s$. According to Lemma 1, we know that the bilinear term $n_{sr}{z}_{sr}$ can be replaced with $n_{sr}$. Then, the constraints can be expressed as follows: ${\sum }_{p\in P\left(r\right)}{w}_{srp}\le Q_s{n}_{sr}$. Let the total amount of CO₂ generated by all emission sources on route $r$ be denoted as follows: $E_r={\sum }_{p\in P\left(r\right)}{E}_p$. From Constraint (6), we know that: ${\sum }_{s\in S}{\sum }_{r:p\in P\left(r\right)}{w}_{srp}+e_p=E_p$, $\forall p\in P$, which ensure that all CO₂ generated by each source is either transported or emitted. Therefore, for each route $r$ which is assigned to ship $s$ ($n_{sr}>0$), the value $E_r$ equals the sum of the CO₂ volume transported by ships from emission sources $p\in P\left(r\right)$ to all storage sites connected to these sources, plus the total CO₂ volume directly emitted from $p\in P\left(r\right)$, that is,

(16)$\sum _{s\in S}\sum _{r\in R}\sum _{p\in P\left(r\right)}{w}_{srp}+\sum _{p\in P\left(r\right)}{e}_p=E_r.$

Since each emission source is included by at least one route, and is served by at least one ship on each route, the following inequality holds:

(17)$\sum _{p\in P\left(r\right)}{w}_{srp}\le \sum _{s\in S}\sum _{r\in R}\sum _{p\in P\left(r\right)}{w}_{srp}.$

Combined with (16), we have:

(18)$\sum _{p\in P\left(r\right)}{w}_{srp}+\sum _{p\in P\left(r\right)}{e}_p\le E_r.$

According to Constraint (11), $e_p$ is non-negative for all $p\in P$, and thus we can conclude that: ${\sum }_{p\in P\left(r\right)}{e}_p\ge 0$. Then, we have ${\sum }_{p\in P\left(r\right)}{w}_{srp}\le E_r$. Therefore, we can replace $Q_s{n}_{sr}$ in the original Constraint (3) with $\mathrm{m}\mathrm{i}\mathrm{n}\left(Q_s{n}_{sr}, E_r\right)$. The modified constraints become:

(19)$\sum _{p\in P\left(r\right)}{w}_{srp}\le \mathrm{m}\mathrm{i}\mathrm{n}\left(Q_s{n}_{sr}, E_r\right), \forall s\in S, \forall r\in R,$

which completes the proof. □

Based on Lemmas 1–4, we establish the following main theorem, which summarizes the equivalence between the original MIP model and its modified formulation incorporating integrality constraints, linearization, and tightened capacity bounds, while preserving the submodular property of the objective function.

Theorem 1.

The original and modified MIP models are equivalent, and the submodular objective function is preserved under this modification.

Let P be the original MIP problem with objective function $F_{\Omega }(\mathit{x},\mathit{z},\mathit{n},\mathit{e})$ and feasible region Ω. Let [P] be the modified MIP problem with objective function ${ F\prime }_{\Omega }(\mathit{x},\mathit{z},\mathit{n},\mathit{e})$ and feasible region Ω’ (defined by (1′)–(11)). Then, [P] is equivalent to [R]. Specifically, the following applies: (i) Any feasible solution to [P] can be transformed into a feasible solution to [R] with the same objective value. (ii) Any feasible solution to [R] can be transformed into a feasible solution to [P] with the same objective value. (iii) Consequently, an optimal solution to [P] corresponds to an optimal solution to [R], and they share the same optimal objective value.

The reformulated model is presented as follows:

[R] min	${\displaystyle \sum _{i\in I}}{C}_i{x}_i+{\displaystyle \sum _{s\in \mathrm{S}}}{\displaystyle \sum _{r\in R}}{C}_s^{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{p}}{z}_{sr}+{\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r\in R}}T\left(F_{sr}+B_{sr}\right)n_{sr}+{\displaystyle \sum _{p\in P}}T\widehat{C}{e}_p$		(20)
s.t.	${\displaystyle \sum _{r\in R}}{z}_{sr}\le 1$	$\forall s\in S$	(21)
	${\displaystyle \sum _{p\in P\left(r\right)}}{w}_{srp}\le \mathrm{m}\mathrm{i}\mathrm{n}\left(Q_s{n}_{sr}, E_r\right)$	$\forall s\in S, \forall r\in R$	(22)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:i\left(r\right)=i}}{\displaystyle \sum _{p\in P\left(r\right)}}{w}_{srp}\le \mathrm{m}\mathrm{i}\mathrm{n}\left(V_i,E_i\right)$	$\forall i\in I$	(23)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:p\in P\left(r\right)}}{n}_{sr}\ge N_p$	$\forall p\in P$	(24)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:p\in P\left(r\right)}}{w}_{srp}+e_p=E_p$	$\forall p\in P$	(25)
	$0\le n_{sr}\le N_{sr}{z}_{sr}$	$\forall s\in S, \forall r\in R$	(26)
	$x_i\in \left\{0,1\right\}$	$\forall i\in I$	(27)
	$z_{sr}\in \left\{0,1\right\}$	$\forall s\in S,\forall r\in R$	(28)
	$w_{srp}\ge 0$	$\forall s\in S, \forall r\in R, \forall p\in P$	(29)
	$0\le e_p\le E_p^{\mathrm{m}\mathrm{a}\mathrm{x}}$	$\forall p\in P$	(30)

4. Experiments

This section conducts computational experiments to verify the effectiveness of our proposed model. Specifically, we implement and solve the MIP model [R], referring to Equations (1)–(11). The experiments were conducted on a desktop computer equipped with 3.40 GHz of 13th Gen Intel Core i7 CPU and 32 GB of RAM, and the MIP model was solved by the Gurobi Optimizer 10.0.1 via the Python 3.11.5 API. In the Gurobi configuration for this model, the MIP gap parameter was explicitly set to 0.1% to specify a relative optimality gap, allowing the solver to terminate when the best feasible solution is within 1% of the optimal bound and balancing solution quality with computational efficiency. Other parameters remained at their default settings, including the utilized number of CPU cores, pre-solve, solver-integrated cuts, and heuristics, and log information. We first set initial values for parameters to obtain basic results. Furthermore, sensitivity analyses were conducted to examine the impact of these parameters.

4.1. Data Collection

The planning horizon is set to 20 years ($T=20$), with annual calculations and constraints applied consistently to facilitate comprehensive long-term planning and system evaluation. A 600 by 600 (n mile) simulation environment is developed, representing a scale pertinent to establishing a maritime CO₂ transport network and analogous to scenarios connecting dispersed industrial clusters to offshore storage opportunities, such as those considered in European CCS planning involving continental sources and potential North Sea storage locations. Within this environment, ten emission sources and five candidate storage sites are uniformly distributed to reflect a typical regional dispersion of major industrial facilities and potential geological storage formations. Distances between any two locations are measured using the Euclidean metric, with $d_{ip}$ denoting the distance between storage site $i$ and emission source $p$, and $d_{p_1{p}_2}$ indicating the distance between two emission sources $p_1$ and $p_2$. For a route including one emission source and one storage site, the route’s sailing distance $l_r$ is calculated as $l_r=d_{ip}+d_{pi}=2d_{ip}$. For a route including two emission sources and one storage site, the route’s sailing distance is calculated as $l_r=d_{i{p}_1}+d_{p_1{p}_2}+d_{p_{2}i}$.

The parameters are categorized into three groups:

Parameters of storage sites. Each candidate storage site is characterized by its annual storage capacity $V_i$, ranging from 10 to 30 million tons, derived from reported capacities of operational CCS projects across diverse geological contexts, such as those assessed in European regions including Denmark and the North Sea [49,50]. The construction cost of storage site $i$ is formulated as $C_i=a_i+b_i{V}_i$. The parameter values are informed by global CO₂ storage initiatives. Specifically, the fixed cost $a_i$ ranges from USD 200 to 600 million, reflecting the investment required for infrastructure development in regions such as the North Sea and Norwegian continental shelf, as documented in global cost assessments and Norwegian CCS projects [50,51]. The variable cost $b_i$ varies between 15 and 40 USD/ton, aligning with reported costs for CO₂ injection into saline aquifers and depleted oil and gas fields, particularly from Norwegian demonstration projects [51]. Specific parameter values for the five candidate storage sites are summarized in Table 4.

Parameters of emission sources. Each emission source has an annual CO₂ emission volume $E_p$, ranging from 2.0 to 5.0 million tons, consistent with reported ranges for major industrial emitters such as coal-fired power plants, cement production facilities, and steel manufacturing units in heavy industrial regions globally, including Europe, the U.S., and China [52,53,54]. The minimum annual berthing requirement $N_p$ is positively correlated with each source’s emission volume. The specific values are detailed in Table 5.

All CO₂ generated by each source is either transported or emitted. The penalty cost for direct atmospheric CO₂ emissions $\widehat{C}$ is set as 110 USD/ton, based on the European Union emission regulation guidelines [55], using a conversion rate of EUR 1 = USD 1.1. This standard applies uniformly across all EU ETS-covered installations, including coastal factories and power plants, as confirmed by the European Commission and national implementing agencies [56,57]. Each source is subject to an upper annual direct emission limit that cannot exceed 5% of its annual CO₂ emission amount, i.e., $E_p^{\mathrm{m}\mathrm{a}\mathrm{x}}=0.05 E_p$. This specific threshold, limiting direct emissions from each CO₂ source to 5% of its total annual CO₂ generation, is established to reflect stringent environmental performance expectations and aligns with ambitious objectives, such as those indicated by the International Maritime Organization’s 2023 GHG Strategy aiming for the adoption of near-zero emission technologies. Furthermore, such high capture rates are supported by lifecycle assessments which emphasize the need for substantial emission reductions to ensure the overall environmental integrity and climate benefit of CCS operations [58,59,60].

Parameters of ships. A fleet of 40 ships is considered, consisting of 10 small ships, 20 medium ships, and 10 large ships. Based on [61], small ships have transport capacities ranging from 2000 to 5000 tons and incur a total chartering cost of USD 20 million over the 20-year planning horizon. Medium ships carry 8000 to 12,000 tons with a 20-year chartering cost of USD 60 million. Large ships transport 18,000 to 25,000 tons and require a total chartering cost of USD 120 million for the same period. Each berthing operation, whether at a storage site or an emission source, incurs a cost of USD 0.05 million for small ships, USD 0.1 million for medium ships, and USD 0.2 million for large ships, in line with 2025 typical port charges for vessels of varying sizes, adjusted for the specialized requirements of CO₂ transport, and supported by BEIS port fee estimates [62]. All ships use very-low-sulfur fuel oil (VLSFO). Referring to [63], the hourly consumption (kg/h) of ship $s$ sailing at speed $v_s$ is expressed as follows:

(31)$f_s\left(v_s\right)=C_{\mathrm{s}}^0+C_{\mathrm{s}}^1·{v_s}^{n_s},$

where $n_s$ is the class-specific speed exponent, which is 3.5 for small ships, 4.0 for medium ships, and 4.5 for large ships. $C_{\mathrm{s}}^0$ and $C_{\mathrm{s}}^1$ are coefficients determined by the tonnage of ship $s$. The parameters for ship classes are presented in Table 6, with $C_s^0$ and $C_s^1$ derived from the median of the value ranges provided in [63]. All ships are assumed to operate at a standardized sailing speed of $v_s=20$ knots. Each year provides $t_0=8500$ h effective operational hours after accounting for downtime and maintenance, allowing for a realistic estimation of the maximum number of annual round trips as $N_{sr}=t_0/\left(l_r/v_s\right)$. The fuel cost for ship $s$ completing one full trip on route $r$ is calculated as $F_{sr}=C_f·l_r/v_s·f_s\left(v_s\right)$, where $C_f$ is the unit fuel price, set as 0.587 USD/kg based on global bunker prices [64].

After establishing the parameter settings, we used these values to derive the basic results, and then we conducted sensitivity analysis to examine the impacts of these parameters.

4.2. Result Interpretation

This section presents the detailed basic results derived from the proposed MIP model, covering storage site construction strategy, capacity utilization of each constructed sites, ship chartering and route assignment strategy, CO₂ transport and emission outcomes, transport operations, and a comprehensive cost analysis.

4.2.1. Storage Site Construction Strategy and Capacity Utilization

In this section, we present the results of the MIP model solved by Gurobi. The results related to the construction strategy are presented in Table 7.

Table 7 summarizes the selected storage sites for construction and their capacity utilization. Among the five candidate sites, $i_1$ and $i_4$ are selected for construction. Site $i_1$, with an annual capacity of 10 million tons, receives 8.9 million tons of CO₂ annually, yielding an 89.0% utilization rate. Site $i_4$ demonstrates even higher efficiency, with an annual intake of 24.7 million tons against a capacity of 25 million tons, achieving a utilization rate of 98.8%. The average utilization rate of the two storage stations is 93.9%, indicating that the constructed sites efficiently match transportation demands.

4.2.2. Ship Route Assignment and Transport Operations

Based on the optimal layout of CO₂ storage sites derived from the MIP model, an analysis is conducted on the ship chartering, routing assignment, and how frequently each chartered ship completes round trips along its designated route per year.

Table 8 presents detailed results regarding ship chartering and route assignments. A total of 12 ships were selected for chartering from the available fleet, consisting of three small-class, seven medium-class, and two large-class ships. Each ship is assigned to a fixed route, connecting one or two emission sources to a designated storage site. Smaller ships are associated with higher values of $n_{sr}$, indicating more frequent annual operations. This is consistent with their lower transport capacity, which necessitates a greater number of trips to complete the required CO₂ transport. In contrast, larger ships are assigned lower values of $n_{sr}$, as their greater capacity allows them to meet transport demand with fewer trips. Additionally, the large ship $s_{40}$ is assigned to a route that connects two emission sources. Its high capacity enables it to carry CO₂ from two sources in a single round trip, making it particularly suitable for joint service. This strategy enhances ship utilization and reduces the need for additional ships, further contributing to cost reduction.

4.2.3. CO₂ Transport and Emission Analysis

Based on the optimized assignment of ships and transport plan, Table 9 presents the detailed allocation of transported CO₂ volumes and direct emissions for each emission source. Specifically, for each emission source, the table lists the amount of CO₂ transported annually by each assigned ship.

As shown in Table 9, 33.5914 million tons of CO₂ are transported annually, with only 0.0086 million tons emitted directly, corresponding to an average emission rate of less than 0.03%. Notably, all CO₂ in emission sources $p_6$ and $p_9$ are transported without any direct emissions. The extremely low level of direct emissions can be attributed to the relatively high penalty cost, which guides the system toward sustainable emission-reducing decisions.

4.2.4. Cost Components Analysis

Based on the results derived from the MIP model, Table 10 presents a comprehensive breakdown of the total system cost in the planning horizon.

The total cost amounts to USD 17,830 million, comprising expenditures related to storage site construction, ship chartering, fuel consumption, berthing operations, and penalties for direct CO₂ emissions. Among these components, the berthing cost accounts for the largest proportion of the total, reaching USD 14,354 million (80.5%). This high cost is primarily driven by the large number of berthing operations required. The construction cost totals USD 1600 million (9.0%), representing the capital investment needed to ensure adequate capacity of selected storage sites. The fuel cost is USD 1137 million (6.4%), incurred during sailing operations across the network. The chartering cost amounts to USD 720 million (4.0%), corresponding to the long-term lease of twelve ships assigned according to transport demand. The penalty cost for direct CO₂ emissions is only USD 19 million (0.1%), which reflects the model’s strong preference for transporting nearly all CO₂, driven by the relatively high penalty imposed on untransported emissions.

4.3. Sensitivity Analysis

In this section, we further analyze the sensitivity of our MIP model to the changes in input parameters, including the annual available sailing time, the sailing speed of all ships, and the capacities of the candidate storage sites.

The parameters of the experiments are set as shown in Table 11. Generally, we conduct 30 experiments with an experiment ID (EID) indexed from 0 to 29. Each of the experiments is included in the corresponding group with a group ID (GID). G0, G1, and G2 are the groups of experiments that aim to illustrate the performance of the MIP model in solving instances with different annual available sailing time ($t_0$), the sailing speed of ships ($v_s$), and the ratio of the actual capacity of each storage station to its original capacity, respectively. The ratio is denoted as $\alpha$, where the original capacity is $V_i$, the actual capacity is $V_i^{\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$, and $\alpha =V_i^{\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}/V_i$. In the “$t_0$ (h)”, “$v_s$ (knots)”, and “$\alpha$” columns, $[a, b, c]$ represents a list of numbers generated from $a$ to $b$ with a step size of $c$.

4.3.1. Different Annual Available Sailing Time

To investigate the impact of varying the annual available sailing time on the model’s effectiveness, we design instances by altering the annual available sailing time from 3000 h to 8500 h while keeping the sailing speed ($v_s=$ 20) and the capacity ratio ($\alpha =1$) unchanged.

Figure 1 shows the optimal results of the MIP model. The x-axis represents the EID, where instances are arranged in increasing order of the annual available sailing time. The y-axis corresponds to the values of indicators, including the objective of an optimal solution (Obj), the number of storage sites constructed (SSN) and the construction cost (SSC) in the planning horizon, the average utilization rate (AUR) of all storage sites, the number of ships chartered (SCN) and the chartering cost (SCC), the amount of CO₂ emitted directly to the atmosphere (CED) in the planning horizon, and the total transportation cost (TTC) in the planning horizon and the average annual number of round trips (TRN) of all ships. Table 12 lists the explanation, calculation, and practical implication of each indicator.

Figure 1a indicates that extending the annual available sailing time from 3000 to 8500 h leads to a steady decline in Obj, before stabilizing around 8000 h. This trend is primarily driven by reductions in SSC and SCC. This observed decrease in Obj with increased sailing time reflects real-world logistics, where greater operational availability of assets typically enhances system efficiency and reduces overall costs. The model’s ability to capture this fundamental relationship supports its practical applicability. A longer sailing time allows each ship to complete more round trips per year, thereby reducing the number of ships needed to meet the transport demand. Therefore, as the annual available sailing time increases, SCN and SCC are reduced, as shown in Figure 1d. This reduction in SCN and SCC is a consequence of improved fleet utilization, a primary objective in capital-intensive transport operations, thereby aligning the model with practical fleet management goals.

As shown in Figure 1b, SSN remains constant at two across all instances due to fixed CO₂ generation volumes and unchanged site capacities. However, SSC shows notable fluctuations, suggesting that although the number of constructed sites remains unchanged, the specific site selection varies. This variation suggests that a longer annual available sailing time enables the model to reconfigure the spatial allocation of selected sites to better align with optimal-cost transport flows. This dynamic site selection, despite a constant number of sites, demonstrates the model’s robustness and sophistication in adapting infrastructure choices to optimize the entire system as operational conditions change, rather than rigidly adhering to a single site configuration. Such adaptability is crucial for real-world strategic planning where operational parameters can vary.

Figure 1c shows that AUR does not exhibit a monotonic trend with respect to sailing time but instead shows distinct fluctuations across different EIDs. It remains around 61% under shorter sailing times (EID 0–5), rises sharply to over 80% at EID 6, drops back at EID 7, and subsequently rises to a peak from EID 9 onward. This pattern suggests that the variation in AUR is not directly driven by sailing time but is primarily influenced by the capacities of the selected storage sites. The model’s ability to reflect how AUR is more closely tied to the specific characteristics of the chosen storage sites, rather than solely with sailing time, which indicates its proficiency in handling complex interactions among different system components, contributing to robust decision support.

Figure 1e illustrates that CED exhibits an overall increasing trend as sailing time increases, which is closely related to the reduction in SCN. After assigning the majority of CO₂ to ships, a small amount of residual CO₂ may remain at certain sources. Whether these residual volumes are transported depends on a cost comparison: only if the transportation cost for the remaining portion is lower than the corresponding emission penalty will the model choose to transport it. Otherwise, direct emission becomes the cost-minimizing option. As the sailing time increases, each ship can complete more round trips, allowing the model to reduce SCN significantly. Although this improves transport efficiency and reduces SCC, having fewer chartered ships makes it more challenging to allocate ships to specific routes, thereby making it harder for the model to find a ship to transport the remaining residual CO₂ on a route at a cost lower than the direct emission, especially for routes with long distances that incur high fuel costs, or that include two emission sources, leading to high berthing costs. As a result, the model tends to opt for direct emissions more, contributing to the overall rising trend in CED, despite some fluctuations due to changing storage site configurations and chartered ships. This observed increase in CED with reduced SCN highlights a critical real-world trade-off where optimizing for lower fleet and chartering costs can lead to increased direct emissions if transporting small, difficult-to-reach CO₂ volumes becomes economically unviable compared to incurring an emission penalty. The model’s capability to explicitly capture and make decisions based on this economic trade-off between transport cost and emission penalty is a strong indicator of its real-world applicability and its robustness in finding a truly minimal-cost solution considering all factors.

Figure 1f shows that TRN increases steadily with sailing time. Given that the average route length across the system does not change substantially, an extended sailing time naturally allows each ship to complete more round trips within a year. As a result, the system-level TRN rises accordingly. In contrast, TTC remains relatively stable, exhibiting only minor fluctuations. This indicates that despite the increase in TRN, the total number of transport operations and the cumulative sailing distance do not grow significantly. Since TTC is jointly influenced by SCN and TRN, the observed stability reflects a compensatory relationship: as the available sailing time increases, SCN decreases while TRN increases. The opposing movement of these two factors offsets each other, leading to limited variation in overall transportation costs. The model’s demonstration of stable TTC, despite significant changes in SCN and TRN, showcases its ability to correctly balance interconnected cost factors. This realistic interplay, where increased efficiency per asset (higher TRN) allows for fewer assets (lower SCN) while maintaining overall transport activity levels, is a testament to the model’s robust cost accounting and its practical relevance for comprehensive system evaluation.

Overall, extending the annual available sailing time reduces both SCC and SSC, thereby lowering Obj. However, reduced SCN results in increased CED due to difficulty in assigning ships on routes with low transportation costs for residual CO₂ volumes. The results highlight that the model minimizes Obj by jointly optimizing SSC, SCC, TTC and CED.

Table 13 shows the optimal results of the model under varying annual available sailing time, ranging from 3000 to 8500 h. As $t_0$ increases, Obj decreases steadily before stabilizing, indicating that a longer sailing time improves overall cost efficiency by enabling each ship to complete more round trips. Most performance indicators change with increasing sailing time. In contrast, SSN remains constant at two across all instances. This suggests that the initial two storage sites are sufficient to meet the transport demand across all levels of annual available sailing time. Attempting to build a third site would result in diminishing returns, where the additional savings in transport and emission costs would no longer offset the increased construction cost. This phenomenon aligns with Lemma 3, which states that the objective function is submodular with respect to the number of constructed storage sites. The submodular property implies that the marginal benefit of adding a storage site decreases as more sites are introduced. In this case, the first two sites contribute significantly to reducing the transportation cost, but beyond a certain point, additional sites would not provide more benefits.

4.3.2. Different Sailing Speed

To investigate the impact of varying the sailing speed of all ships on the model’s effectiveness, we design instances by altering the sailing speed from 10 knots to 30 knots while keeping the annual available sailing time ($t_0=8500$) and the capacity ratio ($\alpha =1$) unchanged. Figure 2 shows the optimal results of the MIP model.

As shown in Figure 2a, Obj exhibits a U-shaped pattern as the sailing speed increases. At low speeds, more ships are needed to meet the transport demand, leading to high SCN and SCC. As speed increases, each ship completes more round trips, reducing SCN and SCC, as shown in Figure 2d. However, beyond a threshold, these values stabilize due to the routing constraint that limits each ship to serving at most two emission sources. This U-shaped Obj curve is characteristic of real-world speed optimization problems across various transport modes. The model demonstrates this practical reality where initial speed increases yield efficiency gains (reduced SCN/SCC), but excessive speed leads to disproportionately higher operational costs (primarily fuel), rendering the system less economical. Capturing this established trade-off is crucial for the model’s credibility and applicability in determining optimal operational strategies.

According to Figure 2b, SSN remains unchanged across all instances. SSC also stays nearly constant, except for a drop at EID 17, which reflects a change in the selected storage site configuration. This change in configuration concurrently affects AUR, which is depicted clearly in Figure 2c, since each storage site differs in capacity.

Furthermore, Figure 2e illustrates a steady increase in CED with rising sailing speeds. As SCN decreases, routing flexibility becomes more limited. With fewer ships available, it becomes more difficult for the model to find a route for transporting the remaining CO₂ at a cost lower than the direct emission penalty. In addition, higher sailing speeds lead to increased fuel costs on the same routes, making it even harder to economically allocate these remaining residual volumes. As a result, the model opts for direct emissions more frequently, leading to an overall rise in CED. The model’s nuanced handling of CED under varying speeds, attributing its increase to both reduced SCN (limiting routing options) and higher fuel costs per trip (making transport of residual CO₂ more expensive), demonstrates a robust understanding of the interconnected factors influencing emission decisions. This ability to weigh multiple cost drivers is vital for practical environmental and economic assessments.

Figure 2f shows that TRN increases with sailing speed since each trip takes less time, allowing ships to complete more round trips within the same annual available sailing time. However, when the sailing speed exceeds a certain threshold, SCN remains constant, and TRN also stabilizes. TTC does not follow a strictly increasing pattern: it decreases slightly at EID 12–16 due to changes in ship chartering, as ships with lower fuel and berthing costs are preferred. At EID 17–22, however, according to Equation (20), the high exponent values (3.5, 4.0, and 4.5 for the different ship classes in this study) cause disproportionately large increases in fuel consumption, resulting in a clear upward trend in TTC. The model’s sensitivity to the exponential relationship between speed and fuel consumption contributes to its real-world applicability as fuel costs are a dominant factor in maritime operations. The clear upward trend in TTC at higher speeds, directly linked to this exponential fuel use, realistically reflects the severe economic penalties of excessive speed, thereby validating the model’s cost calculations and its ability to identify a practically optimal speed range.

Overall, increasing sailing speed reduces SCN and SCC by allowing each ship to complete more round trips within a year, which helps lower Obj. However, as speed continues to rise, fuel consumption rises significantly, causing TTC to increase and ultimately driving Obj back up. In addition, the reduction in SCN leads to higher CED due to limited routing flexibility. The results highlight that the model minimizes Obj by jointly optimizing SSC, SCC, TTC and CED under varying sailing speeds.

Table 14 presents the optimal results under different sailing speeds from 10 to 30 knots. Obj first decreases and then increases as speed rises, with accompanying changes in SCN, SCC, CED, TRN, and TTC. A shift in SSC and AUR reflects changes in storage site selection. SSN remains constant at two across all instances, which suggests that constructing two storage sites is sufficient to meet the transport demand and achieve significant total cost savings.

4.3.3. Different Capacity Ratio of Storage Sites

To investigate the impact of varying the capacity ratio of all storage sites on the model’s effectiveness, we design instances by altering the capacity ratio from 0.4 to 1.6 while keeping the annual available sailing time ($t_0=8500$) and the sailing speed ($v_s=$ 20) unchanged. Figure 3 shows the optimal results of the MIP model.

As shown in Figure 3a, Obj exhibits an overall decreasing trend as the capacity ratio increases. This is primarily due to the reduction in SSC. With a larger capacity ratio, fewer sites are needed to satisfy the CO₂ storage requirement. As a result, both SSN and SSC decrease as the capacity increases, and eventually only one storage site is selected, as reflected in Figure 3b. However, since higher-capacity sites are more expensive to construct, there are less pronounced decreases or even slight increases in SSC in some cases. The model’s strategic response of reducing SSN as individual site capacities increase is a logical adaptation to changing infrastructure parameters. This ability to make high-level decisions on the scale and number of facilities based on their individual characteristics demonstrates its utility for long-term investment planning in real-world CCS projects. Furthermore, the nuanced behavior where SSC might not always decrease monotonically (if larger sites have a disproportionately higher unit construction cost) indicates the model balances economies of scale against potential cost premiums for larger infrastructure, a common real-world consideration.

Figure 3c shows that AUR fluctuates considerably across EIDs, reflecting differences in the number, location, and capacity of selected storage sites. These factors significantly influence the level of storage utilization under each scenario.

As shown in Figure 3d,f, SCN, SCC, TRN, and TTC remain relatively stable. Given that the sailing time, sailing speed, and annual CO₂ generation are fixed, transport operations do not vary significantly. The slight fluctuations in TTC are mainly driven by changes in routing distance under different storage site configurations. The observed stability in transport-related costs (SCN, SCC, TRN, TTC) when only site capacity is varied (while operational parameters like sailing time and speed are held constant) effectively demonstrates the model’s capability to isolate and analyze the impact of specific infrastructure decisions. This feature significantly enhances the capability for comprehensive scenario analysis in practical planning, enabling decision-makers to systematically evaluate the individual impacts of infrastructure versus operational changes.

Figure 3e highlights noticeable variations in CED, primarily due to differences in storage site construction strategies and routing assignments. Although ship chartering decisions vary only slightly, they still influence the routing flexibility and transport capacity, thereby affecting CED. The model’s indication that CED varies based on storage site strategies, even with stable overall transport operations, emphasizes the interconnectedness of infrastructure choices and their downstream environmental consequences via routing adjustments. This highlights the model’s comprehensive view, which is crucial for integrated system planning.

Overall, when varying the capacity ratio, the variation in Obj is mainly driven by changes in SSC. Other cost components remain relatively stable, and CED variation has a limited impact on total cost. These results demonstrate the model’s ability to balance construction, chartering, transport, and emission costs to achieve the minimum system cost.

Table 15 shows the numerical results under varying capacity ratios. Obj gradually decreases as the capacity increases, primarily due to the reduction in SSC, as SSN drops from four to one. AUR varies across EIDs, consistent with differences in the number and location of selected sites. Other indicators, including SCN, SCC, TRN, and TTC, remain relatively stable due to fixed sailing time and speed. CED exhibits relatively noticeable variation, reflecting the influence of storage site and routing strategies, although its absolute values remain small across all cases.

To sum up, the numerical experiments, encompassing variations in annual available sailing time, sailing speed, and storage site capacity ratio (as detailed in Table 13, Table 14 and Table 15), demonstrate that the model not only effectively responds to changes in key parameters by adjusting storage site construction, ship chartering, and route assignment strategies, but also does so in a manner that reflects real-world operational logic and economic trade-offs. The results reveal how different operational and infrastructure configurations influence system-wide outcomes, particularly the Obj and its major components—SSC, SCC, and TTC. Furthermore, these analyses highlight the model’s robustness in navigating complex interactions, such as balancing fleet efficiency against emission targets, or optimizing sailing speed against escalating fuel costs. The model’s ability to identify these critical thresholds and trade-offs, leading to holistic, minimal-cost solutions that consider strategic infrastructure investments, operational efficiencies, and environmental implications simultaneously, supports its practical applicability as a decision-support tool for designing and managing sustainable maritime CO₂ transport systems.

The computational efficiency of the proposed MIP model was assessed by recording the CPU time required to solve each scenario in the sensitivity analyses. As summarized in Table 13, Table 14 and Table 15, the average solution times across the three sets of sensitivity analyses were 71.4 s, 58.9 s, and 39.5 s, respectively. These findings suggest that the model is computationally tractable and can be effectively solved within reasonable time, thereby confirming its practical computational efficiency.

5. Two-Stage Optimization Framework and Solution Approach

In practical scenarios, employing the comprehensive MIP model from Section 3 Problem Formulation can present significant challenges. This is primarily because the model, with its numerous binary and continuous variables and complex operational constraints governing site selection, ship assignment, and routing, inherently belongs to the class of NP-hard optimization problems. For such problems, the computational effort to find a guaranteed optimal solution can grow exponentially with network scale, such as the number of sources, sites, ships, and routes. This inherent complexity means that for large-scale, real-world instances, achieving exact optimal solutions with the integrated MIP model alone can become computationally prohibitive within practical timeframes, despite its comprehensive nature. Moreover, applying such a complex model often requires considerable mathematical optimization expertise. Therefore, to address these issues, particularly the potential computational intractability for larger systems due to its NP-hard nature and the need for a more accessible decision framework, a simplified two-stage modeling approach is proposed. This approach still addresses CO₂ storage site construction, ship chartering, route assignment, and transport operation planning by decomposing the problem into strategic site selection and tactical ship assignment and routing stages. The goal is to provide more tractable and adaptable solutions, especially for real-world applications where scalability and ease of use are critical.

The reasons for employing a two-stage approach are as follows:

Firstly, significantly reduced computational time for large-scale instances: While an MIP model provides optimality, its solution time can become excessively long for complex, real-world scenarios characteristic of NP-hard problems. Heuristic approaches, like the GA used in our second stage, are designed to find high-quality solutions much more rapidly, making them viable for larger problems or when quick decision-making is paramount.

Secondly, enhanced accessibility and reduced reliance on expensive commercial solvers: Our two-stage approach is structured such that the first stage (site selection) can potentially be addressed using simpler, perhaps even greedy, methods to obtain good initial configurations, or other less computationally intensive optimization techniques. The second stage (ship assignment and routing) then employs a genetic algorithm. Crucially, this entire two-stage framework can be implemented and executed without necessitating a commercial MIP solver like Gurobi. This is a significant practical advantage for many organizations, particularly smaller companies or those in regions with limited access to such software, for whom the substantial cost of commercial solver licenses (which can be tens or even hundreds of thousands of dollars) is a prohibitive barrier. Open-source solvers, while available, often exhibit significantly lower computational efficiency for complex MIPs compared to their commercial counterparts. Our two-stage heuristic approach therefore offers a cost-effective and more broadly accessible pathway to sophisticated decision support.

5.1. Two-Stage Modeling Approach

In this section, we describe the structure of the proposed two-stage modeling approach, which aims to simplify the original MIP model by decoupling strategic infrastructure decisions from operational transport planning. The first stage determines the optimal selection of CO₂ storage sites to ensure sufficient total capacity, while the second stage focuses on optimizing ship routing and transport operations given the selected sites.

5.1.1. First Stage: Strategic Site Construction

The first stage is formulated as a binary optimization model. It identifies which candidate storage sites should be constructed to minimize total construction costs while ensuring that the selected sites together offer enough capacity to store all generated CO₂, as shown below:

min	${\sum }_{i\in I}{C}_i{x}_i$		(32)
s.t.	${\sum }_{p\in P}{E}_p\le {\sum }_{i\in I}{V}_i{x}_i$		(33)
	$x_i\in \left\{0,1\right\}$	$\forall i\in I$	(34)

The objective function (32) minimizes the construction cost of selected CO₂ storage sites. Constraint (33) ensures the total amount of CO₂ generated annually by all emission sources does not exceed the total capacity of all selected storage sites. Constraint (34) is the domain of decision variables.

The storage site construction decisions determined in this stage are denoted by the set $I*$, representing the subset of candidate locations $I$ that are selected for construction, i.e., $I*\subseteq I$, $x_i=1$, $\forall i\in I*$. Notably, only the storage sites $i\in I*$ will be considered in the second stage.

Given the parameters of emission sources and candidate storage sites, the first-stage model offers a straightforward yet effective framework for identifying cost-efficient CO₂ storage site combinations. As a binary optimization problem with limited decision variables, it is computationally lightweighted and can be solved rapidly, even in large-scale scenarios. Moreover, the decision logic—selecting a subset of sites to satisfy total storage demand while minimizing cost—is highly intuitive and easily interpretable by industry practitioners. This makes the model particularly well suited for use in early-stage planning, feasibility studies, or stakeholder discussions where transparency and ease of implementation are critical. The selected sites from this stage form the foundation for the second-stage transport assignment model, which determines the ship assignments and routing strategies accordingly.

5.1.2. Second Stage: Operational Ship Assignment and Route Planning

Given the selected storage sites, this stage optimizes the chartered ships, assignment of routes, number of annual round trips, and transport volumes, with the objective of minimizing total tactical costs including chartering cost, transportation cost, and direct emission penalties. The model in this stage can be written as follows:

min	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r\in R}}{C}_s^{ship}{z}_{sr}+{\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r\in R}}T\left(F_{sr}+B_{sr}\right)n_{sr}+{\displaystyle \sum _{p\in P}}T\widehat{C}{e}_p$		(35)
s.t.	${\displaystyle \sum _{r\in R}}{z}_{sr}\le 1$	$\forall s\in S$	(36)
	${\displaystyle \sum _{p\in P\left(r\right)}}{w}_{srp}\le \mathrm{m}\mathrm{i}\mathrm{n}\left(Q_s{n}_{sr}, E_r\right)$	$\forall s\in S, \forall r\in R$	(37)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:i\left(r\right)=i}}{\displaystyle \sum _{p\in P\left(r\right)}}{w}_{srp}\le \mathrm{m}\mathrm{i}\mathrm{n}\left(V_i,E_i\right)$	$\forall i\in I^{*}$	(38)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:p\in P\left(r\right)}}{n}_{sr}\ge N_p$	$\forall p\in P$	(39)
	${\displaystyle \sum _{s\in S}}{\displaystyle \sum _{r:p\in P\left(r\right)}}{w}_{srp}+e_p=E_p$	$\forall p\in P$	(40)
	${0\le n}_{sr}\le N_{sr}{z}_{sr}$	$\forall s\in S, \forall r\in R$	(41)
	$z_{sr}\in \left\{0,1\right\}$	$\forall s\in S,\forall r\in R$	(42)
	$w_{srp}\ge 0$	$\forall s\in S, \forall r\in R, \forall p\in P$	(43)
	$0\le e_p\le E_p^{\mathrm{m}\mathrm{a}\mathrm{x}}$	$\forall p\in P$	(44)

The objective function (35) minimizes the total cost, which consists of three components: the total chartering cost of all chartered ships, the transportation cost in the planning horizon (including fuel and berthing costs), and the penalty cost for direct CO₂ emissions into the atmosphere. Constraint (36) ensures that each ship can be assigned to at most one route, reflecting the condition that not all ships must be chartered. Constraint (37) denotes that, for each ship assigned to a route, the total amount of CO₂ transported per year does not exceed the product of its transport capacity and the number of annual round trips. Constraint (38) requires that the total injected CO₂ into each constructed storage site per year does not exceed its capacity. Constraint (39) requires that each emission source is visited by ships a sufficient number of times per year to meet its minimum berthing requirement, thereby ensuring an adequate transport frequency and enabling timely removal of the generated CO₂. Constraint (40) ensures that all CO₂ generated by each source is either transported or emitted. Constraint (41) imposes operational bounds on the annual number of round trips that each chartered ship can make on its assigned route, ensuring that this value does not exceed the maximum allowable trips based on ship speed and route distance. Constraints (42)–(44) are the domains of decision variables. Constraints (36)–(44) in this model can be viewed as a reduced and reformulated subset of the constraints originally presented in the integrated MIP model, specifically corresponding to Constraints (2)–(7) and (9)–(11). Notably, since the site selection decision is fixed in the first stage, the original storage capacity constraint (4) is reduced to (38), which only enforces capacity compliance at the predetermined sites. This reformulation leads to a more compact and computationally tractable model, facilitating faster solution times and improved scalability. Moreover, the reduced model structure enhances practical interpretability, making it more suitable for real-world applications where decision-makers require operational guidance under established infrastructure conditions.

5.2. Solution Approach and Experiment Results

In this section, we introduce a heuristic solution method to address the second stage of the problem and present the corresponding computational results for CO₂ transport assignment. Given the selected storage sites from the first stage, we adopt a GA to solve the CO₂ transport assignment problem. The GA provides a metaheuristic-based alternative to exact methods, enabling near-optimal solutions under complex constraints.

5.2.1. Genetic Algorithm Design and Parameter Settings for Second Stage Model

In this section, we develop a heuristic algorithm to solve the second-stage CO₂ transportation planning problem, which focuses on assigning ships to feasible routes and scheduling annual round trips based on the storage site construction outcomes from the first stage. This problem exhibits a complex structure characterized by binary decisions for ship chartering and ship-route assignments, continuous variables for round-trip frequencies, and additional continuous decisions related to the direct emission of untransported CO₂. These decisions are subject to multiple constraints, including storage capacity limits, ship transport capacity limits, ship operation boundaries, and direct emission limits. As the problem size scales with the number of ships and routes, the exact solution through an MIP model is computationally intractable. Therefore, we adopt a metaheuristic approach to efficiently explore the solution space while still generating feasible and interpretable results.

To this end, we design a GA tailored to the characteristics of the problem. In the GA framework, each individual (chromosome) encodes a complete assignment of ships to routes and their corresponding round-trip frequencies, ensuring that only the routes connecting to selected storage sites from the first stage are considered. Each gene in the chromosome represents a specific ship and consists of three key components: the assigned route (comprising one storage site and one or more emission sources), and the number of annual round trips $n_{sr}$. The assigned route also indicates whether the ship is chartered. If no route is specified, the ship is considered not chartered. The transported volume is computed based on the ship’s capacity and its trip frequency and directly influences the amount of untransported CO₂, which is represented as direct emission $e_p$ at each source. This gene structure enables the algorithm to jointly optimize ship chartering, routing, and scheduling decisions while satisfying constraints on transport capacity, storage capacity, source berthing frequency, and CO₂ emission volumes.

To initialize the population, a hybrid seeding approach is adopted to generate both feasible and cost-efficient individuals. First, for each emission source, the algorithm constructs a candidate route pool by ranking all feasible single-source and dual-source routes to each selected storage site in ascending order of expected cost, which includes fuel cost, berthing cost, and penalties for direct CO₂ emission based on maximum ship frequency. Then, the ships are sorted by transport capacity and greedily assigned to the lowest-cost routes that still have available storage capacity and unserved emission demand. For each assignment, the number of round trips is determined by the lesser of the ship’s annual upper bound and the amount required to meet the remaining transport volume on that route. After the initial assignment, a repair and reallocation process is applied to resolve infeasibilities. If any emission source still has unmet transport demand, the algorithm either activates an unchartered ship with adequate capacity or reassigns a ship from an over-served source to the under-served one. Reassignment is accepted only if it leads to a reduction in total tactical cost, accounting for fuel and berthing cost increases, changes in emission penalties, and possible improvements in storage utilization. A feasibility checker is applied to ensure that all individuals comply with constraints on source berthing frequency, round trip number limits, ship transport capacity, and storage site capacity. If any constraint is violated, targeted adjustments such as reducing the round-trip frequency or switching to a nearby low-cost alternative route are applied until the individual becomes feasible. To maintain diversity, 50% of the population is generated by perturbing the greedy solution—shuffling ship orders, adjusting round trips, or replacing routes with alternatives of comparable cost. The remaining 50% are randomly generated from the space of feasible ship–route–trip combinations and subsequently repaired using the same logic. This strategy ensures that all initial individuals are feasible, cost-aware, and structurally diverse, which promotes effective exploration and fast convergence in the early generations. Figure 4 presents a sample chromosome layout where each column represents one candidate ship. The rows indicate the assigned storage site, the associated emission source, the number of round trips, and the corresponding transport volume, respectively. Ships that are not chartered are indicated with ‘NA’ entries in the respective rows.

The fitness function directly evaluates the total system cost, including chartering cost, fuel and berthing cost, and penalty cost for direct emission. The algorithm starts from a randomly generated population of feasible solutions and iteratively improves them through selection, crossover, and mutation. Tournament selection is applied to favor high-quality solutions while maintaining population diversity. Crossover is performed using a two-point method to exchange route assignments between chromosomes. Mutation is conditionally applied after the crossover step, specifically targeting individuals ranking in the lowest 30% of the offspring’s fitness values to enhance exploration in poorly performing regions of the solution space. For each selected individual, one ship gene is randomly chosen and mutated with a fixed probability of 0.08. The mutation involves either replacing the assigned route with another feasible alternative from the candidate route pool or adjusting the annual round-trip frequency within its allowable bounds. If the mutation causes any constraint violation, a constraint-based local repair strategy is triggered to restore feasibility without altering the rest of the chromosome. These mutated individuals remain within the same generation and are directly evaluated along with other offspring. The two-point crossover and mutation operations are illustrated in Figure 5 and Figure 6, respectively.

The offspring population replaces the current one, and each chromosome’s fitness is evaluated. The best individuals are carried forward into the next generation, and the cycle of selection, crossover, mutation, and evaluation repeats until either the maximum number of iterations is reached or the best fitness value remains unchanged for a predefined number of consecutive generations.

5.2.2. Computational Results and Analysis

In this section, we present the results of applying GA to solve the second-stage CO₂ transport assignment problem. Given the parameters in Table 5, the total annual CO₂ generation from all emission sources is 33.6 million tons. Among all feasible combinations, selecting storage sites $i_2$ (15 million ton capacity) and $i_4$ (25 million ton capacity) provides a total capacity of 40 million tons, which satisfies the demand. This combination yields the lowest total construction cost of USD 1525 million. This result serves as the foundation for the second stage, where ship assignments are planned to connect all emission sources with the selected storage sites $i_2$ and $i_4$. To solve this assignment problem, the GA is terminated after 300 generations, or earlier if the best fitness value shows no improvement over 50 consecutive generations. The algorithm employs a population size of 100, a crossover rate of 0.8, and a mutation rate of 0.1. Elitism is incorporated by retaining the top 20% of solutions in each generation. These parameter settings are selected to strike a balance between exploration and exploitation and demonstrated to be effective in producing high-quality solutions within an acceptable computational time. The results related to storage capacity utilization are presented in Table 16.

Table 16 shows the two storage sites that were selected during the first-stage construction decision, along with their annual CO₂ receipts and corresponding capacity utilization rates. Site $i_2$, with an annual capacity of 15 million tons, receives 10.6 million tons of CO₂ annually, yielding an 70.7% utilization rate. Site $i_4$ demonstrates even higher efficiency, with an annual intake of 23.0 million tons against a capacity of 25 million tons, achieving a utilization rate of 92.0%. The average utilization rate of the two storage stations is 81.4%, indicating that the GA solution is capable of generating assignments that effectively utilize the capacity of selected storage sites.

Based on the storage site layout determined by the first-stage decision, an analysis is conducted on the ship chartering, routing assignment, and how frequently each chartered ship completes round trips along its designated route per year.

Table 17 presents detailed results regarding ship chartering and route assignments. A total of 12 ships were selected for chartering from the available fleet, consisting of two small-class, seven medium-class, and three large-class ships. Each ship is assigned to a fixed route, connecting one or two emission sources to a designated storage site. Smaller ships are associated with higher values of $n_{sr}$ and larger ships are assigned lower values of $n_{sr}$, as their greater capacity allows them to meet transport demand with fewer trips. Additionally, the large ships $s_{38}$, $s_{39}$, and $s_{40}$ are assigned to routes that connect two emission sources. Compared with the result shown in Table 8, the GA solution charters one more large ship and one fewer small ship. This adjustment may be attributed to the change in storage site selection, which leads to longer average route distances and thus necessitates larger transport capacity per trip to maintain operational efficiency.

Based on the optimized assignment of ships and transport plan by solving the second-stage using GA, Table 18 presents the transport volumes and direct emissions for each emission source.

As shown in Table 15, 33.5873 million tons of CO₂ are transported annually, with only 0.0127 million tons emitted directly, corresponding to an average emission rate of less than 0.04%. Notably, all CO₂ in emission source $p_6$ are transported without any direct emissions. Compared to the results in Table 7, the direct emission volume under the GA solution is slightly higher. This may be attributed to the absence of one small ship with a lower transportation cost and the fact that the average transport route lengths under the GA assignment are generally longer than those in our previously proposed MIP model.

Based on the results obtained through the two-stage optimization framework and the assignment determined by solving the second stage using GA. Table 19 presents a comprehensive breakdown of the total system cost in the planning horizon.

Compared with the cost structure presented in Table 10, the results detailed in Table 19 show no significant change. The berthing cost remains the dominant contributor, accounting for over 80% of the total cost in both solutions. However, slight increases are observed in the chartering cost, fuel cost, and berthing cost under the GA-based assignment. The two-stage approach employing GA resulted in a total cost of USD 18,266 million, which is 2.4% higher than that achieved using the exact MIP model. In terms of computational efficiency, the GA completed the second-stage assignment in an average of 28.3 s, representing a 57.9% reduction compared to the 67.2 s required by the Gurobi-solved MIP model for equivalent operational planning. These results indicate that, for the tested problem scale, the customized GA offers a notable relative decrease in computational time for the operational subproblem with only a marginal increase in cost. While the MIP model proved effective for this particular instance, the improved speed of the GA, along with its independence from costly commercial solvers and reduced demand for specialized optimization expertise, makes it a more practical option. This is especially true in scenarios where rapid decision-making is critical, the problem size increases beyond the capacity of exact methods, or resource constraints limit the deployment of sophisticated MIP tools.

6. Conclusions

Maritime CO₂ transport plays an important role in supporting CCS systems by connecting emission sources with suitable CO₂ storage sites. However, the transport process faces several planning and operational challenges, particularly when long sailing distances and strict emission constraints are considered simultaneously. The locations of CO₂ storage sites significantly affect total tactical costs, as improper site selection may lead to increased sailing distances and higher transportation costs. Once storage sites are determined, further challenges arise in transport assignment, including the selection of appropriate ship types, the allocation of ship-route assignments, and the scheduling of annual sailing operations for each ship. Therefore, optimizing both the siting of CO₂ storage facilities and the assignment of transport tasks is essential to reduce total tactical costs and ensure the effective operation of maritime CCS systems.

To address these challenges, we propose an MIP model for the CSSL-TA problem. We establish a series of lemmas to analyze the mathematical structure of the model. Specifically, we prove that the decision variables associated with CO₂ storage site construction and ship-route assignments must be binary since fractional construction or half-chartered ships are not practical. We also prove that bilinear terms involving ship activity can be replaced by linear expressions, reducing model complexity. Furthermore, we demonstrate that the objective function is submodular with respect to construction, chartering, and routing variables, which reflects the diminishing marginal benefits of adding more resources. Additionally, we tighten two key sets of capacity constraints by coupling supply and demand. Instead of using fixed capacity parameters, the tightened constraints dynamically reflect the actual CO₂ generation volumes associated with each storage site and transport route. These enhancements help reduce the feasible space and accelerate convergence. These theoretical insights provide a solid foundation for enhancing solution efficiency and extending the model in future research.

In the numerical experiments, the model yields a total tactical cost of USD 17,830 million in the planning horizon. Two storage sites are selected for construction, with an average utilization rate of 93.9%. Out of a fleet of 40 available ships, 12 ships are selected for chartering, including three small-class, seven medium-class, and two large-class ships. Among all cost components, the berthing cost accounts for the largest share, reaching 80.5% of the total, while the penalty cost for direct CO₂ emissions is the smallest at only 0.1%. This outcome results from the model’s preference to avoid the relatively high penalty imposed on untransported emissions. The strategy not only ensures sufficient storage capacity and effective transport coverage, but also minimizes costs by aligning storage sitting and ship assignment decisions with capacity constraints, operational needs and direct emission limits. The sensitivity analysis reveals that variations in annual available sailing time, sailing speed, and storage capacity affect the total cost, storage site deployment, chartering cost, transportation cost, and direct emission volume. A longer sailing time, higher sailing speed, and greater capacity lead to fewer ships being needed. These factors, in turn, affect storage site construction strategies, transportation costs, and the volume of direct CO₂ emissions to varying degrees. In particular, higher sailing speeds cause a rapid increase in fuel cost, which significantly raises the total cost and limits the model’s ability to assign ships to routes with low transportation cost for residual CO₂.

To enhance practical applicability, we further design a two-stage planning framework that separates strategic site selection and transport assignment. In the first stage, storage sites are selected based on storage capacity and construction cost. In the second stage, a GA is applied to solve the transport assignment problem given the fixed site locations. While the two-stage approach provides a computationally efficient alternative, the resulting total tactical cost is 2.4% higher than that obtained by the exact MIP model, but it reduces the computational time by 57.9%. This result demonstrates that the proposed two-stage method can serve as a reliable and efficient decision-support tool, particularly useful for real-world applications where decision-makers must balance accuracy and tractability.

The optimization framework aligns with a multi-level decision-making structure:

Strategic level: The framework supports long-term decisions on the construction and spatial deployment of CO₂ storage sites. These decisions involve large-scale infrastructure planning and capital investment, and are typically made by government agencies, regional CCS authorities, or public infrastructure coordinators.

Routing and assignment level: The framework determines which transport routes should be activated, assigns specific ships to selected routes, and specifies the number of round trips for each ship. These mid-level planning tasks are usually managed by CCS system operators or logistics planners responsible for transport design and resource allocation.

Operational level: It provides detailed support for day-to-day execution, including CO₂ volume allocation from emission sources, compliance with berthing frequency requirements, and control of untransported emissions. These operational tasks are implemented by maritime transport companies or CCS logistics service providers.

While the proposed optimization framework demonstrates strong potential for improving the cost effectiveness and feasibility of maritime CO₂ transport, its real-world implementation may be affected by several practical risks. These include uncertainties in CO₂ emission volumes and storage capacities, unplanned infrastructure failures at storage sites or ports, and disruptions in ship operations due to weather conditions, scheduling conflicts, or route congestion. To mitigate such risks, the framework supports proactive measures such as the use of redundant storage capacity, flexible route selection, adaptive ship assignments, and compliance-based emission control strategies. In the event of unforeseen failures, the model allows for real-time operational adjustments and the use of backup transport resources to minimize the severity of consequences. Future research should aim to more explicitly incorporate the impact of realistic factors such as weather conditions, route congestion, and berthing delays, which are crucial in maritime CO₂ transport and can significantly impact the feasibility and reliability of optimization models in real-world applications. The development of stochastic programming or robust optimization models would allow for a more direct representation of these uncertainties, leading to solutions that better ensure operational feasibility and reliability under variable conditions. In addition, the current model can be extended by incorporating real-world complexities such as uncertainty in CO₂ capture volumes, fuel prices, and system performance, which are critical to the robustness of maritime CCS planning. Recent studies have highlighted how stochastic optimization, Monte Carlo simulations, and scenario-based models can capture uncertainties in injection rates, storage efficiency, and energy market volatility, enabling more adaptive infrastructure design and avoiding stranded CO₂ scenarios [65,66,67]. Another promising extension involves integrating flexible pipeline transport and marine shipping options for hybrid CO₂ routing. Research comparing techno-energetic characteristics of liquefied CO₂ shipping and pipelines has demonstrated that such configurations can improve system resilience and reduce costs in maritime contexts [68]. Additionally, the model could incorporate policy-driven constraints such as dynamic emission limits or carbon-trading schemes. Emerging studies in maritime transport operations under emission-trading frameworks show that fluctuating carbon prices significantly influence routing, operational strategies, and investment planning [69]. Lastly, the model can be enhanced by adopting multi-objective optimization approaches that jointly consider cost efficiency, emission reduction, and infrastructure utilization. In particular, formulations that incorporate uncertainty can support more balanced decision-making and help identify Pareto-efficient infrastructure strategies [70]. These extensions would further enhance the model’s applicability to long-term sustainable planning of maritime CCS systems.

Author Contributions

Conceptualization, S.W.; methodology, Y.T., Y.Y., Y.D. and S.W.; software, Y.T.; validation, Y.T.; formal analysis, Y.T.; investigation, Y.T., Y.Y. and S.W.; resources, S.W. and Y.D.; data curation, Y.T. and Y.Y.; writing—original draft preparation, Y.T., Y.Y., S.W. and Y.D.; writing—review and editing, Y.T., Y.Y., S.W. and Y.D.; visualization, Y.T.; supervision, S.W.; project administration, S.W.; funding acquisition, S.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Optimal results of the MIP model with different annual available sailing time (EID: 0–11).

Figure 2 Optimal results of the MIP model with different sailing speeds (EID: 12–22).

Figure 3 Optimal results of the MIP model with different capacity ratios (EID: 23–29).

Figure 4 Structure of a GA chromosome.

Figure 5 Two-point crossover operation (Note: Purple cells represent the original values; yellow cells represent the updated values after crossover operation.).

Figure 6 Mutation operation (Note: Purple cells represent the original values; yellow cells represent the updated values after mutation operation.).

Table 1

Key studies on CO₂ storage site selection.

Paper	Considered Factor						Applicable Scope
	Reservoir Depth	Porosity and Permeability	Caprock Integrity	Storage Capacity	Injectivity Potential	Containment Risk
Bachu [31]	✓	✓	✓	✓			Deep saline aquifers, depleted hydrocarbon fields
Grataloup et al. [32]	✓	✓	✓	✓			Paris Basin
Alcalde et al. [33]			✓	✓	✓	✓	General
Raza et al. [34]	✓	✓		✓	✓		General
Hsu et al. [35]			✓			✓	South Yellow Sea Basin
Luo et al. [36]	✓	✓		✓			General
Aarnes et al. [37]		✓		✓			Commercial-scale storage
Akhurst et al. [38]		✓		✓			European

Note: (1) The “✓” in the table indicates that the respective paper considers the corresponding objective costs; (2) Injectivity Potential: the ability of a reservoir to accept and accommodate CO₂ injection over time, influenced by factors like permeability, porosity, and reservoir pressure; (3) Containment Risk: the risk of CO₂ leakage from the storage site, primarily due to inadequate caprock integrity or potential migration pathways, which could compromise long-term storage security; (4) General: the study’s methods or criteria can be applied broadly to various geological settings (e.g., saline aquifers or depleted fields) without being tied to a specific region or project.

Table 2

Key studies on maritime CO₂ transport (✓ indicates the study considered the corresponding objective or modeling method).

Paper	Objective				Modeling Method
	CO2 Transport Cost	CO2 Direct Emission Cost	Ship Chartering Cost	Storage Site Construction Cost
Zhang et al. [9]	✓			✓	MIP
Roussanaly et al. [12]	✓		✓		Techno-economic modeling (TEM)
Elahi et al. [40]	✓			✓	MIP
Knoope et al. [41]	✓			✓	Real Options Analysis (ROA)
Gabrielli et al. [42]	✓			✓	MIP
Engebø et al. [43]	✓				Fault tree analysis, Monte Carlo simulations
Kjärstad et al. [44]	✓			✓	TEM
Morbee et al. [45]	✓				Node clustering
Bennæs et al. [46]	✓		✓		MIP
Our Work	✓	✓	✓	✓	MIP, two-stage optimization used GA

Table 3

Notations used in the model’s formulation.

Sets
$I$	The set of candidate locations where CO2 storage sites can be constructed
$P$	The set of CO2 emission sources
$S$	The set of ships transporting CO2
$R$	The set of feasible CO2 transport routes
Indices
$i$	Index for locations in $I$
$p$	Index for CO2 emission sources in $P$
$s$	Index for ships in $S$
$r$	Index for routes in $R$
Parameters
$T$	The number of years in the planning horizon
$V_i$	The capacity of the $i$-th CO2 storage site
$C_i$	The cost of constructing a candidate CO2 storage site at location $i$
$E_p$	The amount of CO2 produced by source $p$ per year
$E_p^{\mathrm{m}\mathrm{a}\mathrm{x}}$	The maximum amount of CO2 that emission source $p$ is allowed to release into the atmosphere per year
$N_p$	The minimum annual berthing requirement at source $p$
$C_s^{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{p}}$	The chartering cost of ship $s$
$Q_s$	The transport capacity of ship $\mathrm{s}$
$F_{sr}$	The fuel cost for ship $\mathrm{s}$ completing one full trip on route $r$
$B_{sr}$	The berthing cost for ship $\mathrm{s}$ on route $r$
$N_{sr}$	The maximum number of round trips that ship $s$ can make on route $r$ per year
$\widehat{C}$	The penalty cost of emitting a ton of CO2 into the atmosphere
Decision Variables
$x_i$	Binary variable, which equals 1 if the $i$-th candidate CO2 storage site is constructed, and 0 otherwise
$z_{sr}$	Binary variable, which equals 1 if ship $s$ is chartered and assigned to route $r$, and 0 otherwise
$n_{sr}$	Continuous variable, indicating the number of round trips that ship $s$ makes on route $r$ per year
$w_{srp}$	Continuous variable, indicating the amount of CO2 transported from source $p$ on route $r$ by ship $s$ per year
$e_p$	Continuous variable, indicating the amount of CO2 released into the atmosphere by source $p$ per year

Table 4

The parameters for the five candidate storage sites.

Storage Site	${\mathit{V}}_{\mathit{i}}$ (Million Tons)	${\mathit{a}}_{\mathit{i}}$ (USD Million)	${\mathit{b}}_{\mathit{i}}$ $(\mathbf{USD}/\mathbf{Ton})$	${\mathit{C}}_{\mathit{i}}$ $(\mathbf{USD} \mathbf{Million})$
$i_1$	10	400	40	800
$i_2$	15	200	35	725
$i_3$	20	500	25	1000
$i_4$	25	300	20	800
$i_5$	30	600	15	1050

Table 5

The parameters for the ten emission sources.

Emission Source	${\mathit{E}}_{\mathit{p}}$ (Million Tons)	${\mathit{N}}_{\mathit{p}}$
$p_1$	2.0	5
$p_2$	2.2	5
$p_3$	2.5	6
$p_4$	2.7	6
$p_5$	3.0	7
$p_6$	3.5	8
$p_7$	3.8	9
$p_8$	4.2	10
$p_9$	4.7	11
$p_{10}$	5.0	12

Table 6

The parameters for the three class of ships.

Ship Class	${\mathit{C}}_{\mathit{s}}^0$	${\mathit{C}}_{\mathit{s}}^1$	${\mathit{n}}_{\mathit{s}}$	Transport Capacity (Ton)	${\mathit{C}}_{\mathit{s}}^{\mathbf{s}\mathbf{h}\mathbf{i}\mathbf{p}}$ (USD Million)	${\mathit{b}}_{\mathit{s}}$ (USD Million)	Sailing Speed (Knots)
Small	598.65	0.0198	3.5	2000–5000	20	0.05	20
Medium	649.65	0.0040	4.0	8000–12,000	60	0.10	20
Large	600.45	0.0009	4.5	18,000–25,000	120	0.20	20

Table 7

The storage sites selected for construction and their capacity utilization.

Storage Site	${\mathit{V}}_{\mathit{i}}$ (Million Tons)	The Amount of CO2 Received Annually (Million Tons)	$\mathbf{Capacity} \mathbf{Utilization} (\%)$
$i_1$	10	8.9	89.0
$i_4$	25	24.7	98.8

Table 8

Charter of ships and route assignment.

Chartered Ship	$\mathbf{Ship} \mathbf{Class}$	TransportCapacity (Ton)	Emission SourcesIncluded in the Route	Storage SiteIncluded in the Route	${\mathit{n}}_{\mathit{s}\mathit{r}}$
$s_8$	Small	4200	$p_{10}$	$i_4$	301.4
$s_9$	Small	4500	$p_1$	$i_4$	444.4
$s_{10}$	Small	5000	$p_6$	$i_4$	700.0
$s_{24}$	Medium	10,800	$p_8$	$i_1$	113.0
$s_{25}$	Medium	11,000	$p_2$	$i_4$	200.0
$s_{26}$	Medium	11,200	$p_4$	$i_4$	241.0
$s_{27}$	Medium	11,400	$p_5$	$i_4$	263.0
$s_{28}$	Medium	11,600	$p_8$	$i_1$	256.7
$s_{29}$	Medium	11,800	$p_7$	$i_4$	266.5
$s_{30}$	Medium	12,000	$p_{10}$	$i_4$	311.1
$s_{39}$	Large	24,200	$p_9$	$i_4$	194.2
$s_{40}$	Large	25,000	$p_3$, $p_7$	$i_1$	126.1

Table 9

Transported and directly emitted CO₂ volumes for each emission source.

Emission Source	Chartered Ship	${\mathit{w}}_{\mathit{s}\mathit{r}\mathit{p}}$ $(\mathbf{Million} \mathbf{Tons})$	${\mathit{e}}_{\mathit{p}}$ (Million Tons)
$p_1$	$s=9$	1.9996	0.0004
$p_2$	$s=25$	2.1998	0.0002
$p_3$	$s=40$	2.4991	0.0009
$p_4$	$s=26$	2.6990	0.0010
$p_5$	$s=27$	2.9980	0.0020
$p_6$	$s=10$	3.5000	0
$p_7$	$s=29$	3.1452	0.0014
	$s=40$	0.6534
$p_8$	$s=24$	1.2204	0.0019
	$s=28$	2.9777
$p_9$	$s=39$	4.7000	0
$p_{10}$	$s=8$	1.2659	0.0008
	$s=30$	3.7333
Sum	\	33.5914	0.0086

Table 10

The value and proportion of each component of the cost in the planning horizon.

Cost	Value (USD Million)	Percentage (%)
Total Construction Cost	1600	9.0
Total Chartering Cost	720	4.0
Total Fuel Cost	1137	6.4
Total Berthing Cost	14,354	80.5
Total Penalty Cost for Direct CO2 Emissions	19	0.1
Sum	17,830	100

Table 11

Experiment parameter settings.

GID	EID	${\mathit{t}}_0$ (h)	${\mathit{v}}_{\mathit{s}}$ (Knots)	$\mathit{\alpha }$
G0	0–11	[3000, 8500, 500]	20	1
G1	12–22	8500	[10, 30, 2]	1
G2	23–29	8500	20	[0.4, 1.6, 0.2]

Table 12

Key performance indicators.

Indicator	Explanation	Calculation	Practical Implications
Obj	The total system cost in the planning horizon.	Equation (1): Output of the MIP model’s objective function.	Primary measure of overall economic efficiency and the success of the optimization strategy.
SSN	The total count of CO2 storage sites selected for construction from the candidate locations.	${\sum }_{i\in I}{x}_i$	Indicates the scale of new storage infrastructure required; impacts capital investment and land use.
SSC	The total cost for constructing all the selected storage sites.	${\sum }_{i\in I}{C}_i{x}_i$	Quantifies the capital investment in storage facilities; a major component of system costs.
AUR	Average of annual capacity utilization of all constructed CO2 storage sites.	${\sum }_{i\in I}{\displaystyle \frac{{\sum }_{s\in S}{\sum }_{r:i\left(r\right)=i}{\sum }_{p\in P\left(r\right)}{w}_{srp}}{V_i}} /{\sum }_{i\in I}{x}_i$	Measures the efficiency of capital investment in storage capacity.
SCN	The total count of ships chartered from the available fleet to transport CO2.	${\sum }_{s\in S}{\sum }_{r\in R}{z}_{sr}$	Reflects the size of the operational fleet required; influences overall chartering costs and operational complexity.
SCC	The total cost chartering all selected ships.	${\sum }_{s\in \mathrm{S}}{\sum }_{r\in R}{C}_s^{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{p}}{z}_{sr}$	A major component of system costs.
CED	The total amount of CO2 emitted directly to the atmosphere in the planning horizon.	${\sum }_{p\in P}T{e}_p$	Direct measure of the system’s environmental performance and failure to meet capture targets; incurs financial penalties.
TTC	The total transportation cost (including fuel cost and berthing cost) in the planning horizon.	${\sum }_{s\in S}{\sum }_{r\in R}T\left(F_{sr}+B_{sr}\right)n_{sr}$	Reflects the cost for transporting CO2 from emission sources to storage sites.
TRN	The average number of round trips completed per year by each chartered ship.	${\displaystyle \frac{{\sum }_{s\in S}{\sum }_{r\in R}{n}_{sr}}{{\sum }_{s\in S}{\sum }_{r\in R}{z}_{sr}}}$	Indicates the operational intensity and utilization efficiency of the chartered fleet.

Table 13

Optimal results of the MIP model with different annual available sailing time (EID: 0–11).

EID	${\mathit{t}}_0$ (h)	Obj (USD Million)	SSN	SSC (USD Million)	AUR (%)	SCN	$\mathbf{SCC} (\mathbf{USD} \mathbf{Million})$	CED (Million Tons)	TRN	$\mathbf{TTC} (\mathbf{USD} \mathbf{Million})$	CPU Time (s)
0–11	3000	19,145	2	1850	61.0	22	1440	0.046	88.9	15,850	45.4
	3500	18,849	2	1850	61.3	20	1260	0.068	99.6	15,732	71.1
	4000	18,632	2	1850	61.4	18	1080	0.081	117.9	15,693	57.5
	4500	18,486	2	1850	61.5	17	1020	0.073	123.9	15,608	160.2
	5000	18,320	2	1850	61.7	16	960	0.085	152.3	15,501	95.3
	5500	18,177	2	1850	61.8	14	920	0.126	162.1	15,393	57.9
	6000	18,097	2	1775	82.2	14	880	0.113	180.5	15,430	52.2
	6500	18,027	2	1850	62.0	13	860	0.138	209.2	15,302	52.0
	7000	18,015	2	1775	82.5	13	820	0.127	205.3	15,406	41.8
	7500	17,992	2	1600	93.6	12	800	0.198	219.8	15,570	98.1
	8000	17,843	2	1600	93.7	12	760	0.190	240.9	15,462	58.8
	8500	17,830	2	1600	93.9	12	720	0.172	248.8	15,491	67.1
Avg	5750	18,285	2	1775	73.0	15	960	0.118	170.8	15,537	71.4

Table 14

Optimal results of the MIP model with different sailing speed (EID: 12–22).

EID	${\mathit{v}}_{\mathit{s}}$ (Knots)	Obj (USD Million)	SSN	SSC (USD Million)	AUR (%)	SCN	$\mathbf{SCC} (\mathbf{USD} \mathbf{Million})$	CED (Million Tons)	TRN	$\mathbf{TTC} (\mathbf{USD} \mathbf{Million})$	CPU Time (s)
12–22	10	18,549	2	1850	61.0	17	1060	0.068	102.8	15,632	76.0
	12	18,193	2	1850	61.2	16	920	0.089	124.7	15,413	69.8
	14	17,937	2	1850	61.3	16	820	0.126	163.8	15,253	57.3
	16	17,852	2	1850	61.3	14	780	0.130	182.5	15,208	81.5
	18	17,801	2	1850	61.5	13	720	0.149	219.7	15,215	52.8
	20	17,830	2	1600	93.9	12	720	0.172	248.8	15,491	78.0
	22	17,866	2	1600	93.7	11	640	0.192	253.7	15,605	43.4
	24	18,026	2	1600	93.6	10	640	0.219	260.6	15,762	62.2
	26	18,196	2	1600	93.5	10	640	0.247	268.2	15,929	58.0
	28	18,530	2	1600	93.1	10	640	0.284	262.0	16,259	34.4
	30	18,803	2	1600	93.0	10	640	0.321	265.8	16,528	34.0
Avg	20	18,144	2	1714	78.8	13	747	0.182	213.9	15,663	58.9

Table 15

Optimal results of the MIP model with different capacity ratios (EID: 23–29).

EID	$\mathit{\alpha }$	Obj (USD Million)	SSN	SSC (USD Million)	AUR (%)	SCN	$\mathbf{SCC} (\mathbf{USD} \mathbf{Million})$	CED (Million Tons)	TRN	$\mathbf{TTC} (\mathbf{USD} \mathbf{Million})$	CPU Time (s)
23–29	0.4	18,247	4	2390	93.0	13	620	0.291	233.9	15,205	24.8
	0.6	17,916	3	1985	80.0	12	660	0.127	246.6	15,257	33.7
	0.8	17,895	2	1660	76.6	12	660	0.218	251.4	15,551	37.2
	1.0	17,830	2	1600	93.9	12	720	0.172	248.8	15,491	56.4
	1.2	17,605	1	1140	93.3	13	820	0.230	233.1	15,620	54.6
	1.4	17,376	1	1000	96.0	13	780	0.255	235.9	15,568	34.4
	1.6	17,473	1	1100	84.0	13	780	0.248	236.1	15,566	35.6
Avg	1.0	17,763	2	1554	88.1	13	720	0.220	240.8	15,465	39.5

Table 16

First-stage storage site selection results and associated capacity utilization.

Storage Site	${\mathit{V}}_{\mathit{i}}$ (Million Tons)	The Amount of CO2 Received Annually (Million Tons)	$\mathbf{Capacity} \mathbf{Utilization} (\%)$
$i_2$	15	10.6	70.7
$i_4$	25	23.0	92.0

Table 17

The second-stage charter of ships and route assignment.

Chartered Ship	$\mathbf{Ship} \mathbf{Class}$	$\mathbf{Transport} \mathbf{Capacity} \left(\mathbf{Ton}\right)$	Emission SourcesIncluded in the Route	Storage SiteIncluded in the Route	${\mathit{n}}_{\mathit{s}\mathit{r}}$
$s_9$	Small	4500	$p_2$	$i_2$	488.9
$s_{10}$	Small	5000	$p_6$	$i_4$	327.0
$s_{24}$	Medium	10,800	$p_1$	$i_4$	185.1
$s_{25}$	Medium	11,000	$p_4$	$i_4$	245.4
$s_{26}$	Medium	11,200	$p_9$	$i_2$	201.4
$s_{27}$	Medium	11,400	$p_9$	$i_2$	214.3
$s_{28}$	Medium	11,600	$p_5$	$i_4$	186.7
$s_{29}$	Medium	11,800	$p_7$	$i_4$	211.3
$s_{30}$	Medium	12,000	$p_{10}$	$i_4$	268.1
$s_{38}$	Large	23,500	$p_6,p_{10}$	$i_2$	155.2
$s_{39}$	Large	24,200	$p_5,p_8$	$i_4$	207.8
$s_{40}$	Large	25,000	$p_3$, $p_7$	$i_4$	152.1

Table 18

The second-stage CO₂ transport and direct emissions by each source.

Emission Source	Chartered Ship	${\mathit{w}}_{\mathit{s}\mathit{r}\mathit{p}}$ $(\mathbf{Million} \mathbf{Tons})$	${\mathit{e}}_{\mathit{p}}$ (Million Tons)
$p_1$	$s=$24	1.9992	0.0008
$p_2$	$s=9$	2.1999	0.0001
$p_3$	$s=40$	2.4979	0.0021
$p_4$	$s=25$	2.6992	0.0008
$p_5$	$s=28$	2.1653	0.0016
	$s=39$	0.8331
$p_6$	$s=10$	1.6349	0
	$s=38$	1.8651
$p_7$	$s=29$	2.4933	0.0012
	$s=40$	1.3055
$p_8$	$s=39$	4.1968	0.0032
$p_9$	$s=26$	2.2562	0.0012
	$s=27$	2.4426
$p_{10}$	$s=30$	3.2168	0.0017
	$s=38$	1.7815
Sum	\	33.5873	0.0127

Table 19

The value and proportion of each component of the cost in the planning horizon derived from the two-stage optimization framework.

Cost	Value (USD Million)	Percentage (%)
Total Construction Cost	1525	8.3
Total Chartering Cost	820	4.5
Total Fuel Cost	1243	6.8
Total Berthing Cost	14,650	80.2
Total Penalty Cost for Direct CO2 Emissions	28	0.2
Sum	18,266	100

References

1. Change, P.C. Global Warming of 1.5 C; World Meteorological Organization: Geneva, Switzerland, 2018.

2. Transforming Industry through CCUS. Available online: https://www.iea.org/reports/transforming-industry-through-ccus (accessed on 3 April 2025).

3. Leeson, D.; Mac Dowell, N.; Shah, N.; Petit, C.; Fennell, P.S. A techno-economic analysis and systematic review of carbon capture and storage (CCS) applied to the iron and steel, cement, oil refining and pulp and paper industries, as well as other high purity sources. Int. J. Greenh. Gas Control; 2017; 61, pp. 71-84. [DOI: https://dx.doi.org/10.1016/j.ijggc.2017.03.020]

4. Chamwudhiprecha, N.; Blunt, M.J. CO₂ storage potential in the North Sea. Proceedings of the International Petroleum Technology Conference; Bangkok, Thailand, 7–9 November 2011.

5. Orchard, K.; Hay, M.; Ombudstvedt, I.; Skagestad, R.; Joos, M.; Nysæter, G.; Sjøbris, C.; Jarøy, A.G.; Durusut, E.; Craig, J. The Status and Challenges of CO₂ Shipping Infrastructures. Proceedings of the 15th International Conference on Greenhouse Gas Control Technologies, GHGT-15; Abu Dhabi, United Arab Emirates, 18 March 2021.

6. Kjärstad, J.; Skagestad, R.; Eldrup, N.H.; Johnsson, F. Ship transport—A low cost and low risk CO₂ transport option in the Nordic countries. Int. J. Greenh. Gas Control; 2016; 54, pp. 168-184. [DOI: https://dx.doi.org/10.1016/j.ijggc.2016.08.024]

7. Knoope, M.M.J.; Ramírez, A.; Faaij, A.P.C. Investing in CO₂ transport infrastructure under uncertainty: A comparison between ships and pipelines. Int. J. Greenh. Gas Control; 2015; 41, pp. 174-193. [DOI: https://dx.doi.org/10.1016/j.ijggc.2015.07.013]

8. Roussanaly, S.; Brunsvold, A.L.; Hognes, E.S. Benchmarking of CO₂ transport technologies: Part II–Offshore pipeline and shipping to an offshore site. Int. J. Greenh. Gas Control; 2014; 28, pp. 283-299. [DOI: https://dx.doi.org/10.1016/j.ijggc.2014.06.019]

9. Zhang, S.; Liu, L.; Zhang, L.; Zhuang, Y.; Du, J. An optimization model for carbon capture utilization and storage supply chain: A case study in northeastern China. Appl. Energy; 2018; 231, pp. 194-206. [DOI: https://dx.doi.org/10.1016/j.apenergy.2018.09.129]

10. Bjerketvedt, V.S.; Tomasgard, A.; Roussanaly, S. Optimal design and cost of ship-based CO₂ transport under uncertainties and fluctuations. Int. J. Greenh. Gas Control; 2020; 103, 103190. [DOI: https://dx.doi.org/10.1016/j.ijggc.2020.103190]

11. Drax Global. Why and How Is Carbon Dioxide Transported?. Available online: https://www.drax.com/carbon-capture/why-and-how-is-carbon-dioxide-transported/ (accessed on 17 May 2025).

12. Roussanaly, S.; Jakobsen, J.P.; Hognes, E.H.; Brunsvold, A.L. Benchmarking of CO₂ transport technologies: Part I—Onshore pipeline and shipping between two onshore areas. Int. J. Greenh. Gas Control; 2013; 19, pp. 584-594. [DOI: https://dx.doi.org/10.1016/j.ijggc.2013.05.031]

13. DNV. Background on CO2 Transport by Ship. Available online: https://www.dnv.com/maritime/jip/ceto/background/ (accessed on 17 May 2025).

14. Chemistry Learner. Carbon Dioxide Phase Diagram Explanation. Available online: https://www.chemistrylearner.com/co2-phase-diagram.html (accessed on 17 May 2025).

15. Al Baroudi, H.; Awoyomi, A.; Patchigolla, K.; Jonnalagadda, K.; Anthony, E.J. A review of large-scale CO₂ shipping and marine emissions management for carbon capture, utilisation and storage. Appl. Energy; 2021; 287, 116510. [DOI: https://dx.doi.org/10.1016/j.apenergy.2021.116510]

16. Aramco. Carbon Capture, Utilization, and Storage. Available online: https://www.aramco.com/en/what-we-do/energy-innovation/advancing-energy-solutions/carbon-capture-utilization-and-storage (accessed on 17 May 2025).

17. Kokal, S.L.; Sanni, M.L.; Alhashboul, A.A. Design and implementation of the first CO₂-EOR demonstration project in Saudi Arabia. Proceedings of the SPE Annual Technical Conference and Exhibition; Dubai, United Arab Emirates, 26–28 September 2016; SPE: Richardson, TX, USA, 2016; pp. 1-18.

18. Aramco. Driving Green Initiatives. Available online: https://china.aramco.com/en/news-media/global-news/2020/driving-green-initiatives (accessed on 17 May 2025).

19. Northern Lights. What We Do. Available online: https://norlights.com/what-we-do/ (accessed on 17 May 2025).

20. DNV. Northern Lights Shows the Way to Seaborne CCS Solutions. Available online: https://www.dnv.com/expert-story/maritime-impact/Northern-Lights-shows-the-way-to-seaborne-CCS-solutions/ (accessed on 17 May 2025).

21. Equinor. Northern Lights: The World’s First Cross-Border CO2 Transport and Storage Facility. Available online: https://www.equinor.com/energy/northern-lights (accessed on 17 May 2025).

22. Liotta, G.; Stecca, G.; Kaihara, T. Optimisation of freight flows and sourcing in sustainable production and transportation networks. Int. J. Prod. Econ.; 2015; 164, pp. 351-365. [DOI: https://dx.doi.org/10.1016/j.ijpe.2014.12.016]

23. Nguyen, T.B.; Leonzio, G.; Zondervan, E. Supply chain optimization framework for CO₂ capture, utilization, and storage in Germany. Phys. Sci. Rev.; 2023; 8, pp. 1685-1711. [DOI: https://dx.doi.org/10.1515/psr-2020-0054]

24. Koh, S.L.; Jia, F.; Gong, Y.; Zheng, X.; Dolgui, A. Achieving carbon neutrality via supply chain management: Position paper and editorial for IJPR special issue. Int. J. Prod. Res.; 2023; 61, pp. 6081-6092. [DOI: https://dx.doi.org/10.1080/00207543.2023.2232652]

25. Leonzio, G.; Bogle, D.; Foscolo, P.U.; Zondervan, E. Optimization of CCUS supply chains in the UK: A strategic role for emissions reduction. Chem. Eng. Res. Des.; 2020; 155, pp. 211-228. [DOI: https://dx.doi.org/10.1016/j.cherd.2020.01.002]

26. Lu, H.; Ma, X.; Huang, K.; Fu, L.; Azimi, M. Carbon dioxide transport via pipelines: A systematic review. J. Clean. Prod.; 2020; 266, 121994. [DOI: https://dx.doi.org/10.1016/j.jclepro.2020.121994]

27. Vakili, S.; Manias, P.; Armstrong, L.M.; Turnock, S.; Teagle, D.A. Technical, economic, and environmental assessment of CO₂ ship transport in carbon capture and storage. J. Environ. Manag.; 2025; 373, 123919. [DOI: https://dx.doi.org/10.1016/j.jenvman.2024.123919]

28. Aspelund, A.; Mølnvik, M.J.; De Koeijer, G. Ship transport of CO₂: Technical solutions and analysis of costs, energy utilization, exergy efficiency and CO₂ emissions. Chem. Eng. Res. Des.; 2006; 84, pp. 847-855. [DOI: https://dx.doi.org/10.1205/cherd.5147]

29. d’Amore, F.; Bezzo, F. Economic optimisation of European supply chains for CO₂ capture, transport and sequestration. Int. J. Greenh. Gas Control; 2017; 65, pp. 99-116. [DOI: https://dx.doi.org/10.1016/j.ijggc.2017.08.015]

30. Ma, Z.; Chen, B.; Pawar, R.J. Phase-based design of CO₂ capture, transport, and storage infrastructure via SimCCS 3.0. Sci. Rep.; 2023; 13, 6527.

31. Bachu, S. Sequestration of CO₂ in geological media: Criteria and approach for site selection in response to climate change. Energy Convers. Manag.; 2000; 41, pp. 953-970. [DOI: https://dx.doi.org/10.1016/S0196-8904(99)00149-1]

32. Grataloup, S.; Bonijoly, D.; Brosse, E.; Dreux, R.; Garcia, D.; Hasanov, V.; Lescanne, M.; Renoux, P.; Thoraval, A. A site selection methodology for CO₂ underground storage in deep saline aquifers: Case of the paris basin. Energy Procedia; 2009; 1, pp. 2929-2936. [DOI: https://dx.doi.org/10.1016/j.egypro.2009.02.068]

33. Alcalde, J.; Heinemann, N.; James, A.; Bond, C.E.; Ghanbari, S.; Mackay, E.J.; Hazeldine, R.S.; Faulkner, D.R.; Worden, R.H.; Allen, M.J. A criteria-driven approach to the CO₂ storage site selection of east mey for the acorn project in the north sea. Mar. Pet. Geol.; 2021; 133, 105309. [DOI: https://dx.doi.org/10.1016/j.marpetgeo.2021.105309]

34. Raza, A.; Rezaee, R.; Gholami, R.; Bing, C.H.; Nagarajan, R.; Hamid, M.A. A screening criterion for selection of suitable CO₂ storage sites. J. Nat. Gas Sci. Eng.; 2016; 28, pp. 317-327. [DOI: https://dx.doi.org/10.1016/j.jngse.2015.11.053]

35. Hsu, C.W.; Chen, L.T.; Hu, A.H.; Chang, Y.M. Site selection for carbon dioxide geological storage using analytic network process. Sep. Purif. Technol.; 2012; 94, pp. 146-153. [DOI: https://dx.doi.org/10.1016/j.seppur.2011.08.019]

36. Luo, D.; Yuan, Y.; Chen, J.; Li, Q.; Liang, J.; Zhao, H. Structural and reservoir characteristics of potential carbon dioxide storage sites in the northern south yellow sea basin, offshore eastern China. J. Mar. Sci. Eng.; 2024; 12, 1733. [DOI: https://dx.doi.org/10.3390/jmse12101733]

37. Aarnes, J.E.; Selmer-Olsen, S.; Carpenter, M.E.; Flach, T.A. Towards guidelines for selection, characterization and qualification of sites and projects for geological storage of CO₂. Energy Procedia; 2009; 1, pp. 1735-1742. [DOI: https://dx.doi.org/10.1016/j.egypro.2009.01.227]

38. Akhurst, M.C.; Quinn, M.; Delprat-Jannaud, F.; Bossie-Codreanu, D.; Hannis, S.; Lecomte, J.C.; Frykman, P. Risk assessment-led pre-characterisation of a site for commercial-scale storage of CO2. Proceedings of the Fourth EAGE CO2 Geological Storage Workshop; Stavanger, Norway, 8–10 April 2014.

39. Gruson, J.F.; Serbutoviez, S.; Delprat-Jannaud, F.; Akhurst, M.; Nielsen, C.; Dalhoff, F.; Iacobellis, S. Techno-economic assessment of four CO₂ storage sites. Oil Gas Sci. Technol. IFP Energ. Nouv.; 2015; 70, pp. 753-766. [DOI: https://dx.doi.org/10.2516/ogst/2014057]

40. Elahi, N.; Shah, N.; Korre, A.; Durucan, S. Multi-stage stochastic optimisation of a CO₂ transport and geological storage in the UK. Energy Procedia; 2017; 114, pp. 6514-6525. [DOI: https://dx.doi.org/10.1016/j.egypro.2017.03.1787]

41. Knoope, M.M.J.; Ramírez, A.; Faaij, A.P.C. The influence of uncertainty in the development of a CO₂ infrastructure network. Appl. Energy; 2015; 158, pp. 332-347. [DOI: https://dx.doi.org/10.1016/j.apenergy.2015.08.024]

42. Gabrielli, P.; Campos, J.; Becattini, V.; Mazzotti, M.; Sansavini, G. Optimization and assessment of carbon capture, transport and storage supply chains for industrial sectors: The cost of resilience. Int. J. Greenh. Gas Control; 2022; 121, 103797. [DOI: https://dx.doi.org/10.1016/j.ijggc.2022.103797]

43. Engebø, A.; Ahmed, N.; Garstad, J.J.; Holt, H. Risk assessment and management for CO₂ capture and transport facilities. Energy Procedia; 2013; 37, pp. 2783-2793. [DOI: https://dx.doi.org/10.1016/j.egypro.2013.06.163]

44. Kjärstad, J.; Morbee, J.; Odenberger, M.; Johnsson, F.; Tzimas, E. Modelling large-scale CCS development in Europe linking techno-economic modelling to transport infrastructure. Energy Procedia; 2013; 37, pp. 2941-2948. [DOI: https://dx.doi.org/10.1016/j.egypro.2013.06.180]

45. Morbee, J.; Serpa, J.; Tzimas, E. Optimised deployment of a European CO₂ transport network. Int. J. Greenh. Gas Control; 2012; 7, pp. 48-61. [DOI: https://dx.doi.org/10.1016/j.ijggc.2011.11.011]

46. Bennæs, A.; Skogset, M.; Svorkdal, T.; Fagerholt, K.; Herlicka, L.; Meisel, F.; Rickels, W. Optimization of a ship-based logistics system for carbon capture and storage. Proceedings of the International Conference on Computational Logistics; Cham, Switzerland, 1–3 September 2022; pp. 44-59.

47. Ajayi, T.; Gomes, J.S.; Bera, A. A review of CO₂ storage in geological formations emphasizing modeling, monitoring and capacity estimation approaches. Pet. Sci.; 2019; 16, pp. 1028-1063. [DOI: https://dx.doi.org/10.1007/s12182-019-0340-8]

48. Schrijver, A. Combinatorial Optimization; Springer: Berlin, Germany, 2003.

49. Anthonsen, K.L.; Christensen, N.P. EU Geological CO₂ Storage Summary; Geological Survey of Denmark and Greenland: Copenhagen, Denmark, 2021; pp. 1-28.

50. Smith, E.; Morris, J.; Kheshgi, H.; Teletzke, G.; Herzog, H.; Paltsev, S. The cost of CO₂ transport and storage in global integrated assessment modeling. Int. J. Greenh. Gas Control; 2021; 109, 103367. [DOI: https://dx.doi.org/10.1016/j.ijggc.2021.103367]

51. Gassnova, S.F. Potential for reduced costs for carbon capture, transport and storage value chains (CCS). Nor. Full-Scale CCS Demonstr. Proj.; 2020; 28, pp. 1-61.

52. Duan, H.; Dong, X.; Xie, P.; Chen, S.; Qin, B.; Dong, Z.; Yang, W. Peaking industrial CO₂ emission in a typical heavy industrial region: From multi-industry and multi-energy type perspectives. Int. J. Environ. Res. Public Health; 2022; 19, 7829. [DOI: https://dx.doi.org/10.3390/ijerph19137829]

53. U.S. Environmental Protection Agency. Sources of Greenhouse Gas Emissions: Overview. Available online: https://www.epa.gov/ghgemissions/sources-greenhouse-gas-emissions (accessed on 3 April 2025).

54. China CO2 Emission Accounts 1997–2015. Available online: https://www.nature.com/articles/sdata2017201 (accessed on 3 April 2025).

55. Directive, E.C. Directive of the European Parliament and of the Council on the Promotion of the Use of Energy from Renewable Sources. Proceedings of the European Parliament and Council Conference on Renewable Energy; Brussels, Belgium, 23 April 2008.

56. EU Action on the EU Emissions Trading System (EU ETS). Available online: https://climate.ec.europa.eu/eu-action/eu-emissions-trading-system-eu-ets/about-eu-ets_en (accessed on 15 May 2025).

57. Stationary Operators-EU ETS Compliance. Available online: https://www.epa.ie/our-services/licensing/climate-change/eu-emissions-trading-system-/emissions-trading-system---stationary-installations/stationary-operators---eu-ets-compliance/ (accessed on 15 May 2025).

58. 2023 IMO Strategy on Reduction of GHG Emissions from Ships. Available online: https://www.imo.org/en/OurWork/Environment/Pages/2023-IMO-Strategy-on-Reduction-of-GHG-Emissions-from-Ships.aspx (accessed on 15 May 2025).

59. Now Available: July 2024 Edition of the DNV Class Rules for Ships and Offshore Units. Available online: https://www.dnv.com/news/now-available-july-2024-edition-of-the-dnv-class-rules-for-ships-and-offshore-units/ (accessed on 15 May 2025).

60. Global Status of CCS 2021-Global CCS Institute. Available online: https://www.globalccsinstitute.com/resources/publications-reports-research/global-status-of-ccs-2 (accessed on 15 May 2025).

61. Harper Petersen-Container. Available online: https://www.harperpetersen.com/container (accessed on 3 April 2025).

62. CO2 Shipping Study Final Report for BEIS. Available online: https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/761762/BEIS_Shipping_CO2.pdf (accessed on 3 April 2025).

63. Du, Y.; Chen, Q.; Quan, X.; Long, L.; Fung, R.Y. Berth allocation considering fuel consumption and vessel emissions. Transp. Res. Part E Logist. Transp. Rev.; 2011; 47, pp. 1021-1037. [DOI: https://dx.doi.org/10.1016/j.tre.2011.05.011]

64. ShipBunker. Global 20 Ports Average Bunker Prices. Available online: https://shipandbunker.com/prices/av/global/av-g20-global-20-ports-average (accessed on 3 April 2025).

65. Wang, Z.; Qi, S.; Zheng, B. Uncertainty analysis of the storage efficiency factor for CO₂ saline resource estimation. Energies; 2024; 17, 1297. [DOI: https://dx.doi.org/10.3390/en17061297]

66. Roussanaly, S.; Ouassou, J.A.; Anantharaman, R.; Haaf, M. Impact of uncertainties on the design and cost of CCS from a waste-to-energy plant. Front. Energy Res.; 2020; 8, 17. [DOI: https://dx.doi.org/10.3389/fenrg.2020.00017]

67. Middleton, R.; Yaw, S. The cost of getting CCS wrong: Uncertainty, infrastructure design, and stranded CO₂. Int. J. Greenh. Gas Control; 2018; 70, pp. 1-11. [DOI: https://dx.doi.org/10.1016/j.ijggc.2017.12.011]

68. Poulopoulos, S.G.; Tsiakaras, P.; Voutetakis, S. CO₂ marine transportation from a techno-energetic perspective. J. Mar. Sci. Eng.; 2024; 12, 12.

69. Zhang, Y.; Du, B.; Wu, W.; Liao, T.W. Sustainable maritime transportation operations with emission trading. J. Mar. Sci. Eng.; 2023; 11, 1647.

70. Spek, M.V.D.; Fout, T.E.; García, M.; Kuncheekanna, V.N.; Matuszewski, M.; McCoy, S.; Morgan, J.C.; Nazir, S.M.; Ramírez, A.; Roussanaly, S. . Uncertainty analysis in the techno-economic assessment of CO₂ capture and storage technologies. Int. J. Greenh. Gas Control; 2020; 100, 103113. [DOI: https://dx.doi.org/10.1016/j.ijggc.2020.103113]