I. IntroductionUnder current air traffic management (ATM) operations, one of the primary causes for flight delays and congestion is that the number of flights (demand) often exceeds the supply of the airspace accommodation (capacity). The effort thereby to achieve demand and capacity balancing (DCB) is typically known as air traffic flow management (ATFM), which is regarded as an enabler of ATM efficiency and effectiveness [1]. Since the last few decades, a number of researchers have focused their activity on the development of network ATFM models to minimise the congestion costs incurred from demand and capacity imbalances in both airports and airspace volumes (sectors). However, it is important to highlight that the existing network ATFM models present significant challenges in computational tractability when dealing with large-scale problems. As pointed out in [2], this could be one of the main reasons, besides the fairness concern, why they have not been successfully transitioned into practice.In response to this issue, the computational challenges have been targeted by researchers through applying different numerical computation methods. For the models based on the well-studied Bertsimas Stock-Patterson model [3], a Dantzig-Wolfe decomposition method was proved highly efficient to solve the block-angular formulation [4,5], where the network traffic flow can be decomposed flight by flight. As reported in [6], the effects of this parallel speed-up were further improved by an implementation of the Graphics Processing Units (GPUs), due to the large amount of cores over the Central Processing Units (CPUs). A similar column generation approach was also adopted in [7], where ground holds, airborne delays, reroutes and cancellations are considered as possible control actions, and the approach was then extended to generate recourse strategies in the presence of uncertain capacity constraints.With regards to the network ATFM models that are based on Eulerian formulation, a dual-decomposition method was used, which decomposes the traffic flow path by path, and each flight path was solved (optimised) independently [8]. The dimension of a large-scale cell transmission model therefore depends only on the number of identified flight