2.1. The liquefaction train

A liquefaction train at an LNG facility is comprised of a hot and a cold section. NG from the gas field enters hot section, operating at above ambient temperature. Here, NG is pre-treated to remove acid gas (carbon dioxide and hydrogen sulfide), water and mercury. The processed NG then enters the cold section at temperature T1, pressure P1 and flow rate Q1, respectively, as shown in Figure 1. Temperature T1 depends on the geographical location and can vary from 25 to 30 °C, pressure P1 usually ranges from 50 to 60 bar whereas mass flow rate Q1 (tons per day, T/d) depends on the availability of NG. There are different designs for the cold section. In the C3MR design (Lim et al., Reference Lim, Choi and Moon2012), the cold section pre-cools NG in C3 kettles from T1 to temperature T2 and subsequently to T3 in the main cryogenic heat exchanger (MCHE) using a mixed refrigerant (MR). MR consists of nitrogen (N2), C1, C2, and C3. T2 usually approaches −30 to −27 °C whereas T3 ranges from −150 to −145 °C depending on a variety of factors such as NG composition, MR composition and pressure, and flow rates of NG and MR. The C3 kettles and MCHE are shell-and-tube heat exchanger units with NG flowing on the tube side, C3 in the kettles, and MR in the MCHE, both on the shell side. The duty required to circulate propane and MR to cool NG from T1 to T3 is provided by two compressors. Figure 1 illustrates compressor 1 (C3) and compressor 2 (MR). Cooling NG from T1 to T3 results in the vaporization of C3 and MR; vapor heat is ejected to the atmosphere by air or water cooler before returning back to C3 kettles and MCHE respectively. When upstream pressure P1 is high, the final cooling to T4 = −163 °C occurs in the flash vessel, where NG from MCHE is flash evaporated at pressure P4 (close to the atmospheric pressure). As a result, the flow Q3 from the MCHE is divided into a vapor stream with flow rate Q5, and a liquid stream with flow rate Q4, the latter to storage tanks maintained at atmospheric pressure. The vapor stream is EFG to the FG pool, whereas the liquid stream is LNG for export. The nature of the flash evaporation process is such that Q5 $ \ll $ Q4 with $ \mathrm{Q}3=\mathrm{Q}4+\mathrm{Q}5 $ to retain mass balance; the temperatures and pressures of the EFG and LNG are similar.

Figure 1. Schematic of the cold section of LNG train. The end flash vessel shown in blue produces end flash gas, used as fuel gas for the facility.

2.2. Replicate trains

As noted in Section 1, LNG facilities often contain replicate trains; Figure 2 shows a schematic for two replicate trains studied in this article. Here, EFG from the end flash vessels of each train is sent to the FG pool along with other sources of FG such as BOG and LBOG. The FG pool supplies the FG to the LNG facility. When there is excess FG, the flare valve is opened and the excess FG is flared. To prevent flaring, the typical practice is to reduce EFG production from both trains equally, since trains are notionally replicates by design.

Figure 2. Schematic of two replicate trains, Train 1 and Train 2, feeding EFG to FG pool besides BOG from LNG tank LBOB from tank in the loading vessel (also shown in blue). When the FG pool has excess FG it is released and flared through the flare valve.

In this work, we take advantage of replicate flash vessels at the LNG facility to minimize flare value opening. The presence of replicate components such as compressors, MCHEs, coolers, and C3 kettles at LNG facilities generally can be similarly exploited for operational improvements.

3.1. Exploratory analysis

We consider the operation of flash vessel units U₁, U₂ of replicate trains Tr₁, Tr₂, with EFG mass flow rates Q5₁, Q5₂. Figure 1 motivates the assumption that Q5 for individual units is dependent on (a) the corresponding flow Q3 of NG from the MCHE to the flash vessel, (b) the outlet temperature T3 of NG from the MCHE to the flash vessel, and (c) flash vessel pressure P4. The “manipulated” (or, in statistical terminology, “treatment”) variables Q3, T3, and P4 can be changed independently, thereby influencing Q5. We anticipate that increasing the values of Q3 and T3 will lead to a higher value of Q5. Conversely, a higher P4 will lead to a lower Q5.

We seek to assess fairly whether Q5 from U₁ and U₂ is similar. Ideally, we would conduct a series of experiments on both units, where the values of Q3, T3, and P4 were set at common values across trains, and differences in Q5 were quantified. However, such experiments are impractical economically for trains in continuous operation. Nevertheless, over the course of the normal operation of the trains in time, the set points of Q3, T3, and P4 for the two trains vary, exploring a domain of typical set points for both trains. We can therefore exploit these historical data to quantify differences in Q5. It is, of course, critical that our assessment is fair, in particular, because the domains of Q3, T3, and P4 might be different for the two trains. Since Q5 depends on Q3, T3, and P4, it is essential that the historical data for both trains is filtered such that the treatment variables Q3, T3, and P4 correspond to similar sets of values across the two units. Concisely in mathematical notation, we wish to compare Q5 $ \mid $ (Q3,T3,P4) conditionally across trains, rather than Q5 marginally. The simple filter condition applied takes the form

(1) $$ \mathrm{LL}\le {X}_1^t/{X}_2^t\le \mathrm{UL}\hskip1em \mathrm{for}\hskip0.5em \mathrm{all}\hskip0.5em \mathrm{of}\hskip0.5em X=\mathrm{Q}3,\mathrm{T}3,\mathrm{P}4 $$

where $ {X}^t $ is the value of $ X $ and time $ t $ , for data sampled every 5 minutes for a period of a contiguous calendar year. We emphasize that the filter considered is applied to all of $ X= $ Q3, T3, and P4. Further, LL indicates a common lower limit for the ratio of manipulated variables across trains, set at 0.98 in this work. UL indicates the corresponding common upper limit, set at 1.02. The effect of filtering manipulated variables is illustrated in Figure 3, for data corresponding to the calendar year 2019. Panels of the figure are scatter plots of X₂ on X₁ for X = Q3, T3, P4, and Q5, with green dots indicating data for time points at which the filter conditions in Equation 1 are satisfied, corresponding to ~10% of the unfiltered sample, over all years of available data. Data for all other time points is shown in blue. Of course, filtering yields subsets of operational data for Tr₁ and Tr₂ of equal size. For reasons of commercial confidentiality, note also that all flow Q3 and Q5 presented in this work (e.g., in Figures 3, 4, and accompanying tables in Appendix A) have been normalized using a common factor $ k $ (i.e., Normalized Flow = $ k $ × Observed Flow) such that the maximum Q5 (over all trains and years) in the filtered data is 100 T/d after normalization. No other variables are normalized.

Figure 3. Scatter plots of Q3, T3, P4, and Q5 across trains Tr₁, Tr₂ of operational data sampled at 5-minute intervals for the year 2019. Values for time points satisfying the filter conditions in Equation 1 are shown in green. All other time points are shown in blue. Values of Q3 and Q5 have been normalized using a common factor so that the global maximum value of Q5 is 100 T/d.

Figure 4. Histograms of filtered Q5₁ (blue) and Q5₂ (red) data per annum, for years 2015–2019. Panel titles indicate the number of observations n retained after filtering. Vertical lines and annotated text give mean values of filtered data. Values of Q5 have been normalized using a common factor so that the global maximum value is 100 T/d.

The quality of NG sourced from upstream wells, distributed to the two trains, varies over time. As the proportion of low boiling point components (e.g. C2, C3, and butane, C4) in NG increases, Q5 production reduces in both U₁ and U₂. Moreover, the performance of LNG trains often exhibits seasonal patterns that can influence Q5; filtering (Equation 1) ensures that the comparison of units is not influenced by season and other external variations of the common NG input to liquefaction. Filtering, therefore, allows us to characterize underlying differences in the operational characteristics of the trains, rather than differences in inputs and operating set points.

Since (replicate) trains are optimized through numerical simulations during design, we expect differences in Q5 across trains to be small. We might, therefore, expect that a comparatively long period of historical data might be required to quantify differences in operational characteristics with confidence: in particular, analysis of filtered data from only 1 year can lead to spurious conclusions. Therefore, here, we analyze historical operational data for the 5-year period 2015–2019. The panels of Figure 4 show histograms of filtered Q5 per annum for the years 2015–2019, for train Tr₁ (blue) and Tr₂ (red). The title of each plot shows the year and number n of filtered 5-minute observations. Vertical blue and red lines and annotated text give sample means of filtered Q5₁ and Q5₂, respectively. Figure 4 suggests for each of the 5 years, that Q5 through U₁ is greater than that through U₂. The difference in sample means ranges from 2.5 to 5.5 T/d. Corresponding tables of summary statistics are provided in Appendix A. It is also interesting that the number of observations retained after filtering is considerably higher in 2018 and 2019 than in 2016 in particular, possibly indicating a more consistent setting of operational conditions across trains in more recent years.

The figure for unfiltered data corresponding to Figure 4 is shown in Figure 5. It is notably difficult to see from the figure that there is a material difference between the operating characteristics of trains Tr₁ and Tr₂. This emphasizes the need to consider the conditional behavior of Q5 given its driver variables Q3, T3, and P4.

Figure 5. Histograms of full unfiltered data for Q5₁ (blue) and Q5₂ (red) per annum, for years 2015–2019. Panel titles indicate the number of observations n retained after filtering. Vertical lines and annotated text give mean values of filtered data. Values of Q5 have been normalized using a common factor so that the global maximum value is 100 T/d.

3.2. Statistical testing

The exploratory analysis above suggests that Q5 from vessel U₁ in train Tr₁ is higher than that from U₂ in train Tr₂. We can quantify this using a statistical hypothesis test to assess whether the population mean $ {\overline{\mathrm{Q}5}}_1 $ of Q5 in train Tr₁ is greater than the corresponding population mean $ {\overline{\mathrm{Q}5}}_2 $ in train Tr₂. To perform this one-sided hypothesis test, we set the null hypothesis H₀ that there is no difference between $ {\overline{\mathrm{Q}5}}_1 $ and $ {\overline{\mathrm{Q}5}}_2 $ , and an alternative hypothesis H $ {}_1 $ that $ {\overline{\mathrm{Q}5}}_1 $ $ > $ $ {\overline{\mathrm{Q}5}}_2 $ . Then we calculate whether there is sufficient evidence in the data to reject the null hypothesis in favour of the alternative. Various parametric and non-parametric tests are suggested in the literature (e.g., Marshall and Jonker, Reference Marshall and Jonker2011) for this purpose. The choice of test depends on the nature of the data and the specific question at hand. Here we use the independent two-sample Student’s t-test, calculating test-statistic $ t $ measuring the difference in population means relative to the variability within the groups using sample data. This test assumes that the variances of the two samples are approximately equal. For samples of random variables $ {X}_1 $ and $ {X}_2 $ with common sample size $ n $ , $ t $ is calculated as

(2) $$ t=\left({\overline{X}}_1-{\overline{X}}_2\right)/{s}_d $$

where $ {\overline{X}}_1 $ and $ {\overline{X}}_2 $ are the sample means for Q5₁ and Q5₂ (from Tables A1 and A2 in Appendix A), and $ {s}_d $ is the standard error of the difference in means given by $ {s}_d^2=\left({s}_1^2+{s}_2^2\right)/n $ , where $ {s}_1^2 $ and $ {s}_2^2 $ are corrected sample estimates for the variance of $ {X}_1 $ and $ {X}_2 $ . $ {s}_d $ can also be written as $ {s}_d^2=2{s}_p^2/n $ , where $ {s}_p $ is an estimate for the pooled standard deviation of the samples given by $ {s}_p^2=\left({s}_1^2+{s}_2^2\right)/2 $ . The test statistic $ t $ follows a t-distribution with $ \nu =2\left(n-1\right) $ degrees of freedom (Evans et al., Reference Evans, Hastings and Peacock2000). This probability distribution generalizes the standard normal distribution: both the t-distribution and standard normal distribution have mean zero and exhibit a bell-shaped curve, but the t-distribution has heavier tails controlled by shape parameter $ \nu $ . Typically, the null hypothesis H₀ is rejected at the $ \alpha =0.05 $ level; this occurs when the value of the t-statistic calculated exceeds a critical value $ {t}_{\mathrm{crit},\nu}\left(1-\alpha \right) $ equal to the $ \left(1-\alpha \right)\times 100 $ $ =95 $ %ile of the t-distribution with $ \nu $ degrees of freedom.

(3) $$ \left({\overline{X}}_1-{\overline{X}}_2\right)/{s}_d-{t}_{\mathrm{crit},\nu }(0.95)>0. $$

Multiplying the left-hand side above by $ {s}_d $ gives $ \left({\overline{X}}_1-{\overline{X}}_2\right)-{s}_d\times {t}_{\mathrm{crit},\nu }(0.95) $ , equal to the lower confidence limit LCL for the difference $ {\overline{X}}_1-{\overline{X}}_2 $ in population means. Rejecting H $ {}_0 $ is therefore also equivalent to estimating LCL > 0. For the total sample $ n>100 $ , $ {t}_{\mathrm{crit},2\left(n-1\right)}(0.95)\approx 1.645 $ , the 95%ile of standard normal distribution, to at least two decimal places; for smaller sample sizes, values of $ {t}_{\mathrm{crit},2\left(n-1\right)}(0.95) $ are provided by standard statistical software.

Table 1 shows the results of significance testing for the difference in population mean duty, $ {\overline{\mathrm{Q}5}}_1 $ - $ {\overline{\mathrm{Q}5}}_2 $ , between trains Tr₁ and Tr₂, annually from 2015 to 2019. In percentage terms, $ {\overline{\mathrm{Q}5}}_1 $ exceeds $ {\overline{\mathrm{Q}5}}_2 $ by some 2.8% to 6.4%.

Table 1. Independent two-sample t-test for population mean difference $ {\overline{\mathrm{Q}5}}_1 $ - $ {\overline{\mathrm{Q}5}}_2 $ per annum. Null hypothesis rejected for each year since LCL > 0. Note that the critical value $ {t}_{\mathrm{crit},\nu }(0.95) $ at infinite sample size is adopted as a good approximation, since $ n>1000 $ throughout

When there is evidence that the variance of the two samples is not equal, we can use Welch’s t-test (Welch, Reference Welch1947) as an alternative to the test above. For the current data, using the corresponding Welch test at $ \alpha =0.05 $ , the null hypothesis of equality of $ {\overline{\mathrm{Q}5}}_1 $ and $ {\overline{\mathrm{Q}5}}_2 $ was also rejected for each of the years 2015 to 2019; see Appendix B for details.

3.3. Regression and adjusted regression plots

For each of units U₁ and U₂ on trains Tr₁ and Tr₂, respectively, in turn, we establish linear regression models for Q5 in terms of Q3, T3, and P4 of the form

(4) $$ \mathrm{Q}5=f\left(\mathrm{Q}3,\mathrm{T}3,\mathrm{P}4\right)+\unicode{x025B} $$

for regression function $ f $ , where $ \unicode{x025B} $ is assumed to be a zero-mean Gaussian random variable with unknown standard deviation. Here, we assume that $ f $ takes the linear form

(5) $$ f(\mathrm{Q}3,\mathrm{T}3,\mathrm{P}4)=a+b\ \mathrm{Q}3+c\ \mathrm{T}3+d\ \mathrm{P}4 $$

for parameters $ a $ , $ b $ , $ c $ , and $ d $ to be estimated. Following DuMouchel (Reference DuMouchel1988), we then use adjusted response or adjusted regression plots to quantify the effects of individual treatment variables (more naturally referred to as covariates in a regression context) in regression models for Q5 in terms of Q3, T3, and P4, for each of trains Tr₁ and Tr₂. In essence, these are generalizations of partial residual and augmented partial residual plots (Mallows, Reference Mallows1986), useful for linear regression models with arbitrary power and interaction terms. The fitted regression function $ \hat{f} $ from Equation 5 is

(6) $$ \hat{f}\left(\mathrm{Q}3,\mathrm{T}3,\mathrm{P}4\right)=\hat{a}+\hat{b}\;\mathrm{Q}3+\hat{c}\;\mathrm{T}3+\hat{d}\;\mathrm{P}4 $$

where $ \hat{\bullet} $ represents an estimate. The corresponding residuals from the regression form the set $ {\left\{{r}^i\right\}}_{i=1}^n $ , with

(7) $$ {r}^i=Q{5}^i-\hat{f}\left(\mathrm{Q}{3}^i,\mathrm{T}{3}^i,\mathrm{P}{4}^i\right)\hskip1em \mathrm{for}\hskip0.24em i=1,2,\dots, \mathrm{n} $$

where $ {\left\{\mathrm{Q}{3}^i,\mathrm{T}{3}^i,\mathrm{P}{4}^i\right\}}_{i=1}^{\mathrm{n}} $ is the set of values of Q3, T3, and P4 in the data sample of filtered data for regression model fitting.

Next, adjusted fit functions are calculated for each of Q3, T3, and P4 in turn. For example in the case of Q3, the adjusted fit function is the average value of $ \hat{f} $ , expressed as a function of Q3, over all n observations in the data sample

(8) $$ {g}_{Q3}(q)=\frac{1}{n}\sum \limits_{i=1}^n\hat{f}\left(q,\mathrm{T}{3}^i,\mathrm{P}{4}^i\right). $$

Similar adjusted fit functions can be derived for each covariate in each train in turn. Finally, the set $ {\left\{\tilde{Q}{5}_{Q3}^i\right\}}_{i=1}^n $ of adjusted response values for Q5 with respect to Q3 is calculated using

(9) $$ \tilde{Q}{5}_{Q3}^i={g}_{Q3}\left(\mathrm{Q}{3}^i\right)+{r}^i\hskip1em \mathrm{for}\hskip0.24em i=1,2,\dots, \mathrm{n} $$

where $ {\left\{{r}^i\right\}}_{i=1}^n $ are the residuals from the full regression (Equation 4). Similar sets of adjusted response values can be calculated for response Q5 with respect to each covariate in each train in turn.

Adjusted response values for Q5 with respect to each of Q3, T3, and P4 are shown in Figure 6, for train Tr₁ (blue) and Tr₂ (red). The anticipated directions of the trends of Q5 with covariates are seen in each case. However, despite the trains being nominally replicates, the magnitudes of gradients are larger for train Tr₁ regardless of covariate. Briefly, Q5 is more sensitive to changes in covariates for train Tr₁. To achieve unit reduction in Q5, the reduction in Q3 (and/or T3) needed in Tr₁ is smaller than that needed in Tr₂. This is potentially a valuable handle with which to reduce the need for flaring.

Figure 6. Adjusted response values for Q5 with respect to Q3, T3, and P4 for U₁ (blue circles) and U₂ (red circles). Corresponding adjusted fit functions $ g $ are shown as black lines.

Note that the adjusted regression methodology is applicable generally, regardless of the form of the regression function in Equation 4.

3.4. Implementation of recommendations

Given the findings above, trials were conducted on the liquefaction trains to evaluate the impact on flaring of different manipulations of set-points of manipulated variables on Tr₁ and Tr₂ end flash units U₁ and U₂. In the first period (“Period 1”), each time the flare valve was on the verge of opening, a common reduction of T3 was made for both trains, followed by a common reduction of Q3 if necessary. In the second period (“Period 2”), preferential treatment was given to Tr₁:T3₁ and was reduced first, followed if necessary by Q3₁, T3₂, and Q3₂ if flaring persisted. P4 was not used as a handle during the trial. Results are shown in Figure 7. Panels show the mean flare valve opening in Periods 1 (left) and 2 (right) as a function of the mean T3 (x-axis) and total Q3 (y-axis). The figure indicates a reduction in High and Medium flare valve opening in Period 2 compared with Period 1, resulting in a reduction of up to 45% in flaring-related CO₂ emissions. Polygons (magenta) in each panel show approximate ranges for mean T3 and total Q3 within which the risk of High or Medium flaring is low. The area of the polygon for Period 2 is considerably wider than for Period 1, indicating that the reduction of T3 and Q3 for Tr₁ before those of Tr₂ is advantageous in reducing FG flaring.

Figure 7. Flare valve opening, ranging from High, Medium to Low for Period 1 (left) and Period 2 (right) as a function of a mean of T3 and a sum of Q3 from Tr₁ and Tr₂. In Period 1, simultaneous and equal reductions were made, first for T3 and subsequently if necessary for Q3, for both trains at the point of flare onset. In Period 2, T3₁ and then Q3₁ (if necessary) were reduced first, followed (if necessary) by T3₂ and Q3₂. Polygons show domains of mean T3 and total Q3 corresponding to low risk of High and Medium flare opening.