## 1 Introduction  
The consumption of natural gas can produce 50% less pollution than other fossil fuels [1], and it accounts for 23.7% of global energy consumption [2]. It has also been reported that natural gas will gradually take the place of coal as an important fuel for power generation [3], and many natural gas burners have been designed and modeled [4; 5]. In these applications of industry-scale utilizations, predictive tools like computational fluid dynamics (CFD) represent an effective and economical approach for the burner design and optimization. However, in terms of the turbulent flow, gas-phase chemical reactions, heat transfer, etc., the accurate numerical simulation of natural gas combustion is usually competitive. Many turbulent models, including the standard \(k\)-\(\varepsilon\) model (SKE) [6], the realizable \(k\)-\(\varepsilon\) model (RKE) [7], the shear-stress transport (SST) \(k\)-\(\omega\) model [8], and the Reynolds stress model (RSM) [9; 10; 11]--even the detached eddy simulation (DES) model [12; 13] and the large eddy simulation (LES) model--embedded in commercial computational fluid dynamics CFD software, can address the turbulence flow well [14; 15]. For turbulence-chemistry interaction modeling in computational fluid dynamics (CFD), many models, including the eddy break up (EBU) [16], the eddy dissipation concept (EDC) [17; 18], the mixture fraction probability density functions (PDF) [19; 20; 21; 22], the conditional moment closure (CMC) models [23; 24], and various types of flamelet models [25; 26; 27], have been developed. An extended version of the EBU approach, known as the EDC, has been developed to incorporate detailed chemical kinetics in turbulent flows. This model has shown adequate predictions for premixed, partially premixed, and non-premixed combustion regimes [26]. As one of the most commonly adopted approaches in the modeling of turbulent reacting flows in the context of the steady/unsteady compressible Reynolds-averaged Navier-Stokes equations (RANS/URANS), which sometimes referred to as Reynolds-averaged simulations (RAS), the EDC was initially developed in the 1970s [17; 18]. The EDC was then formulated as a well-established turbulent combustion closure model in the 1990s\(-\)2000s [28; 29; 30]. The EDC has been successfully applied to the numerical simulation of combustion for a long period of time [30; 31], and it is implemented in most available commercial CFD codes. Various models [15; 32; 33; 34; 35; 36; 37; 38; 39] have been validated with the Flame D setup [40; 41; 42; 43] with the detailed chemical kinetic mechanism of GRI-Mech 3.0 [44] (consisting of 53 species and 325 elemental reactions, abbreviated as 53-325) or the reduced mechanisms based of GRI-Mechs 1.2 [45] and 2.11 [46] under steady-state species assumptions. For example, sp21 and sp24 [47; 48] are from GRI-Mech 1.2, and two augmented reduction mechanisms, ARM 9 [49] and ARM 19 [49], are reductions of the GRI-Mech 2.11. In addition to the work of GRI-Mechs, other mechanisms have been also developed for methane-air combustion simulations, such as the skeletal mechanism [50; 51] (6-41), KH97 [52] (28-104), COMB [53] (49-236), K97 [54] (57-353), and SKG04 [55] (72-520). Methane-air combustion was once described by a four-step global reaction mechanism with six species and four reactions [32; 56] for engineering applications. Though this global reaction mechanism was proven to strongly predict CO concentrations, it fails in predicting H\({}_{2}\) and H\({}_{2}\)O, especially in fuel-rich conditions [32; 54]. For applications of industry-scale natural gas combustion furnaces, the computational domains are, however, often much bigger (containing more grid cells) than those in experimental burners used to validate combustion models.