Energy analysis and surrogate modeling for the green methanol production under dynamic operating conditions

Green methanol production, based on intermittent renewable energy sources, requires a more flexible operation mode and close integration with other sections, such as, the electrical grid and electrolysis processes. In this study, methanol synthesis and distillation processes (MSD) for pilot-scale green methanol production (corresponding to 22,236 tons/year) were investigated by dynamic modeling, focusing on energy analysis and dynamic characteristics during load change (LC) operations. The dynamic simulation results with a ramp rate of 50% load / h indicated energy efficiencies of 87.7% (at full-load) and 90.2% (at half-load) for the methanol synthesis process, 86.8% (full-load) and 82.4% (half-load) for the methanol distillation process, and 77.1% (full-load) and 75.4% (half-load) for the MSD process. Relatively small fluctuations were achieved with a ramp time of 1 h for the LC operations. Based on the constructed dynamic model, a surrogate modeling for the MSD process was conducted using the nonlinear autoregressive exogenous model (NARX) model, which exhibited good accuracy with the evaluated performance for the testing data of the root-mean-square error (RMSE) = 3.09 × 10-5, mean absolute error (MAE) = 2.30 × 10-4, and R2 = 1.0. The constructed NARX model can be further integrated with models for other sections of the power-to-methanol process.


Introduction
Global warming due to the greenhouse effect is an urgent issue for human civilization. The reduction of CO 2 emissions from human activities is one of the most important steps for achieving the 1.5 • C and 2 • C temperature targets established by the Paris Agreement [1]; this necessitates the green transition of the energy sector. Increasing shares of intermittent renewable energy sources, such as solar and wind power, will be the primary challenge in balancing the electrical grid. Power-to-x technology, where "x" can represent heat [2], power (e.g., batteries and electric vehicles [3]), gases (e.g., hydrogen and methane [4,5]), or liquid fuels. (e.g., methanol and diesel [6,7]), provides potential energy storage solutions to this issue; using this technology, surplus power can be converted to different forms of energy or energy carriers.
In particular, power-to-methanol (hereinafter refer to as P2M) is a technology with significant potential; green methanol is produced from CO 2 (or CO 2 -rich gas) obtained from renewable carbon sources (e.g., biogas) or captured from manufacturing processes involving heavy CO 2 emissions (e.g., cement plants) and H 2 is produced from electrolysis technologies using renewable electricity [8][9][10] instead of conventional fossil fuels [11][12][13]. This process has attracted growing interest since the concept of the "methanol economy" was proposed by Nobel laureate George A. Olah [14,15].
Green methanol production typically involves the processes of hydrogen production, CO 2 capture or biogas production, storage and purification of the gases, methanol synthesis and distillation. The present study investigated the latter two processes. Heterogeneous catalytic fixed-bed reactors with commercial Cu/Zn/Al 2 O 3 catalysts are widely employed in conventional methanol industries under typical operating conditions of 50-100 bar and 200-300 • C [16,17]. The improvement of these industries is crucial for the methanol synthesis process. The main reactions for this process are described by the following equations: Copper-based commercial catalysts also exhibit good potential for green methanol production based on CO 2 or CO 2 -rich carbon sources [18]. Heterogeneous catalytic methanol synthesis has been employed in project demonstrations, for example, MefCO2 [10] in Germany and Power2Met [18] in Denmark, and a small-scale commercial methanol plant (4000 tons/year) in Iceland [19]. Larger-scale green methanol production has been recently proposed, by the Danish company REintegrate ApS (8000 tons/year) [20], North-C-Methanol project in Belgium (45000 tons/year) [21], and Norwegian companies Statkraft, Finnfjord, and CRI(100000 tons/year) [22].
With the growing deployment of green methanol (MeOH) production, its integration with the electrical grid (also called the P2M process; Fig. 1) based on intermittent renewable energy sources such as solar and wind power brings new challenges. More flexible system operations are required for green methanol production (also for other green chemical production e.g., ammonia [23]), which necessitates more frequent dynamic operations for the methanol synthesis and distillation processes [24][25][26][27]. To overcome the challenges related to the dynamic operations for green methanol production, two primary aspects need to be researched: (1) the dynamic characteristics of the green methanol production system (e.g., energy efficiencies and operational constraints during different dynamic operations) and (2) system integration between green methanol production and the electrical grid (e.g., decisionmaking of plant operations according to the different conditions of the renewable energy sources).
Factors at different scales ranging from catalyst pellets, reactors, and distillation columns to process [28,29], may influence the dynamic characteristics of the complex methanol production system. The present study focuses on the reactor and process aspects through a modeling study. The dynamic modeling of conventional methanol production has not been widely studied. Rezaie et al. [30] conducted long-term (4 y) dynamic simulations of an industrial methanol reactor, focusing on the influence of catalyst deactivation. Abrol  numerical stability of the mathematical model. However, there remains a lack of dynamic modeling for green methanol production at either the reactor or system scale, and the present study aims to fill this gap. Different modeling strategies can be applied for system integration of the different sections shown in Fig. 1. (1) An overall system model can be developed in a single modeling tool based on the first-principle methods on one modeling platform, for example, with Simulink/Matlab; however, this approach is time-intensive because of the complexity of the methanol production system. (2) The sections can be co-simulated using different modeling tools with data communications, for example, cosimulation of Simulink/Matlab and process engineering software such as Aspen Plus. This approach can utilize the advantage of different modeling tools; however, appropriate solutions for simultaneous data exchange are required. (3) The strategy used in this study is based on employing modeling tools such as those in (2) but generating a surrogate model for the methanol production section based on the simulation results by process engineering software, that is, Aspen Plus; subsequently, the surrogate model is combined with the models of other sections on other modeling platforms. This approach presents an alternative method for data exchange between the modeling tools for the sections by using a surrogate model, which can be easily realized by methods such as machine learning and provide a faster computational speed than that based on the first-principle methods [33].
Surrogate models are commonly used to replace complex underlying models and reduce computational effort in chemical and process engineering [33][34][35]. Additionally, they are employed under the condition where there are data regarding input-output relationships but without underlying models [33]. For dynamic processes, McBride and Sundmacher [33] reviewed various applications (e.g., chemical reactions and pharmaceutical processes) of surrogate modeling, involving methods such as kriging and artificial neural networks (ANN). Shokry et al. [35] presented several surrogate models for dynamic chemical processes based on the nonlinear autoregressive exogenous model (NARX), where different metamodel types (ANN and ordinary kriging) were integrated   into the NARX models. Case studies on the dynamic processes involving a bioreactor, three-tank, and oil-shale pyrolysis batch reactor were conducted based on the design of the computer experiment, and the NARX models exhibited good prediction accuracies. Tahkola [36] demonstrated four surrogate models (linear autoregressive model, NARX, Long Short-Term Memory, and Gated Recurrent Unit) for the dynamic power-to-gas process, and the accuracies achieved by all four models seem favorable. Khaleghi et al. [37] developed a NARX model for online health diagnostics of lithium-ion batteries. A high estimation accuracy was achieved by the model, using the health indicators extracted from the voltages as input data.
In the present study, the processes of methanol synthesis and distillation (hereafter referred to as MSD) for green methanol production were investigated by dynamic modeling, which was conducted using the software Aspen Plus Dynamics V9. The load change (LC) operations of the system and key performance indicators (KPIs), such as the energy efficiency and dynamic characteristics of the MSD process, were investigated. Furthermore, surrogate modeling for the MSD process was conducted using the NARX model based on the data generated by the constructed dynamic model.

Dynamic model of the MSD process
The dynamic modeling of the MSD process was conducted using Aspen Plus Dynamics, in which the initial state of the process was extracted from steady-state simulations using the Aspen Plus software. The flowsheet of the MSD process including the methanol synthesis and distillation processes, is shown in Fig. 2.

Methanol synthesis
For the methanol synthesis process shown in Fig. 2, the feed gases of H 2 and CO 2 (provided from other sections) were mixed with a molar ratio of H 2 /CO 2 = 3 (stoichiometric value), and further mixed with the recycle stream (REC). Before being introduced into the methanol reactor (R1), the feed gas (S1) was preheated to 220 • C by heat exchangers HX1 and HX3. HX3 is mainly used for the system start-up (not included in this study), and its duty is zero under the investigated conditions. The hot effluent stream was cooled by heat exchangers HX2 and HX1. The heat exchanged from HX2 was used to heat the reboiler of the distillation column (D1) in the methanol distillation section. The stream (S32) after HX1 was further cooled by CL to 30 • C. Crude methanol (PRO2, containing methanol and water) was obtained in the liquid phase of the gas-liquid separator (SEP and SEP1), which was further purified in the distillation section. Most of the gas phase (VAP) was recycled by a recycle compressor (CMP). A very small ratio (0.1% of VAP) of purging gas (stream PUR1) was set, which aims to minimize the CO 2 emission for the green methanol production. A value lower than 0.1% may cause convergence problem. Using a larger purge ratio can result in lower flow rate of the recycle stream as well as a smaller reactor size but a higher CO 2 loss, e.g., 5.8% of the feed CO 2 is purged with a purge ratio of 1%. The optimal design of this value may need further investigations for both the economic and environmental aspects, which is not included in this study.
A quasi-adiabatic 1D plug-flow model was employed for the fixedbed methanol reactor (R1). The adiabatic fixed-bed reactor showed good potential for green methanol production, which has a lower capital cost than other options, for example, water-cooled and gas-cooled reactors [18]. The heat capacities and weight of the catalyst bed and reactor vessel were considered, which could influence the rate of temperature change due to thermal inertia. The reaction kinetics presented by Vanden Bussche and Froment [38] were employed to calculate the reaction rates under different conditions (validated in Fig. S1 in Appendix A. Supplementary data). The Predictive Soave-Redlich-Kwong (PSRK) model was used as for the equation of state considering the both the methanol synthesis process in the methanol reactor and the twophase system in the gas-liquid separator (SEP) with polar molecules at the high operating pressure, which also showed a better fit when compared to the RKS-MHV2 model [39] and reasonable agreement with the experimental data regarding CO 2 solubility in methanol/CO 2 /water system ( Fig. S2 in Appendix A. Supplementary data). An average effectiveness factor of 0.7 for reactions (2) and (3) was approximately evaluated according to a previous study [18]. A typical catalyst tablet of Φ6 mm × 4 mm (e.g., Topsøe MK-121) was packed in the catalyst bed, and the pressure drop through the catalyst bed was calculated using the Ergun equation [40] (Appendix A. Supplementary data). Table 1 lists the primary operating and reactor parameters for the methanol synthesis section; these parameters are based on the scalingup design (22236 tons/year methanol production) for the pilot plant in the Power2Met project [18]. The hydrogen flow rate for the feed gas is close to 6000 m 3 /h corresponding to electrolysis power consumption of approximately 25 MW. The volume of the main pipeline for this section was considered and added to the design of the SEP volume in the model. A heat transfer coefficient of 0.5 W/m 2 •K was considered for the heat loss of the reactor to the environment, assuming a well-insulated reactor; the heat losses for other parts were not considered. The heat transfer coefficient and heat exchanger area for HX1 were calculated using Aspen EDR (Appendix C. Supplementary data). The heat transfer coefficients for R1 (shown in Table 1) were approximated using correlations (Appendix C. Supplementary data). The coefficients obtained Signal generators (RAMP1-2) were used for dynamic operations, to alter the rates of flow change (i.e., LC) for the feed gases; furthermore, ten proportional-integral (PI) controllers were used for this section, as shown in Fig. 2 and Table 2. The proportional and integral gains were set based on the Ziegler-Nichols and Tyreus-Luyben tuning rules by using the automatic controller tuning in Aspen Plus Dynamics.

Methanol distillation
The crude methanol (stream CRD in Fig. 2) from the methanol synthesis section primarily contained water and methanol with a mole ratio of approximately 1:1, and a small amount of dissolved light gases (mainly CO 2 ). To obtain high-purity methanol, processes with 2-4 distillation columns are commonly employed in conventional methanol plants. For green methanol production, there are much less impurities in the crude methanol [41], and processes with 1-2 distillation columns are sufficient to produce fuel-grade methanol products [19]. This study employed a single distillation column with 28 equilibrium stages. The crude methanol (stream CRD) was preheated to 75 • C by HX4 before entering the distillation column D1. The effluent stream from the top of D1 was cooled by HX5 to 55 • C, where most of the methanol product (MOH) in the gas phase condensed. The light gases were purged after separation from the MOH in the liquid phase in the drum (SEP2). Part of the liquid phase was fed back to the top of the column, and the remaining portion was obtained as the MOH. Waste water (WST) was obtained as the bottom product of the column. The main operating and process parameters for the methanol distillation section are listed in Table 3. Seven PI controllers were used for dynamic operations in this section, as shown in Fig. 2 and Table 4. The purities of the products MOH and WST were mainly adjusted by controlling the flow rate of MOH (i.e., reflux ratio) for the distillation column through valve 13, and the heat duty of the reboiler. The heat duty of the reboiler was adjusted by controlling the temperature (87.6 • C) of the sensitive stage (22). Additionally, lags of 1 min for the temperature measurement were assumed for the control loops [42] in the MSD process.

Surrogate modeling
The framework for developing the surrogate model for the MSD process is shown in Fig. 3, which is based on the dynamic model (step 1) described in Section 2.1. For further integration with other sections in the P2M system, such as electrolysis and electrical grid, the MSD process can be considered as a load for the KPIs, such as energy efficiency and    6 methanol production. The input variables and output KPIs (step 2) are provided in Table 5 (Section 3.2). Steps 3 and 4, shown in Fig. 3, are also described in Section 3.2.
The NARX model is a type of recurrent dynamic neural network [43] that has been widely used for the prediction of time-series data [43][44][45][46]. The NARX model predicts the target output using the following equation [45]: where the output ŷ(t +1) is the predicted output value, F is the nonlinear mapping function, and d u and d y are the time delays for the input u (•) and output y (•) values, respectively. The detailed architecture of the NARX neural network is shown in Fig. 4, where IW, LW, b, and f represent the input weight, layer weight, bias, and activation function, respectively, and TDL is the tapped delay lines of both the exogenous input u (•) and output y (•) values, which are expressed by [46]: The input portion of the NARX model includes the past and present values of u (•) and past value of a real output target y (•) or the predicted target feedback ŷ (•); this generates the two architectures shown in Fig. 5 with open-and closed-loops, respectively. In the present study, the NARX model was developed using an open-loop architecture, including training, validation, and testing. The developed model can be transformed into a closed-loop form for predictions that are multi-step ahead of the target output; furthermore, the closed-loop state can be switched back to the open-loop form when real target output values are known.
The commonly used Levenberg-Marquardt function was selected for the training of the NARX model. The model performance was evaluated by the RMSE, MAE, and R 2 , which are given by the following equations [37]: where N represents the number of data points, and Ŷ i and Y i are the predicted and real target output values, respectively.

Energy analysis
The MSD process is a part of green methanol production (as shown in Fig. 1), which can be considered as a load in the P2M system and evaluated by KPIs such as energy efficiency, upper and lower power limits, and ramp limits. In the following sections, the energy efficiency and other performance indices for the MSD process under different dynamic operations are investigated.
A schematic diagram of the power inputs and outputs for the MSD process is shown in Fig. 6; the hydrogen produced by renewable electricity and methanol produced represent the main power input and output, respectively. Other power inputs include heat exchangers, recycle compressor, and pumps. The purging gases from stream PUR1, which contain a high content of hydrogen, were assumed to be used for combustion and heating other systems. The heat duty of the reboiler was provided by the heat from the HX2 and the possible utilities. The power required for pumping the cooling water for the MSD process was not considered. The energy efficiencies for the methanol synthesis (MS), methanol distillation (MD), and MSD processes were calculated by the following expressions, respectively: η MS =ṁ CRD LHV CRD +ṁ PUR1 LHV PUR1 − P HX2 m H2 LHV H2 + P HX3 + P CMP (6) η MD =ṁ MOH LHV MOḢ m CRD LHV CRD + (P Reb + P HX4 ) + (P P1 + P P2 ) η MSD =ṁ MOH LHV MOH +ṁ PUR1 LHV PUR1 m H2 LHV H2 + (P HX2 + P HX3 + P Reb + P HX4 ) + (P CMP + P P1 + P P2 ) where ṁ (kmol/h) is the mole flow rate of the streams, LHV is the lower heating value for the gases, and P represents the heat duty of the heat exchangers (the duty is negative for cooling process, e.g. HX2) or the power inputs for the recycle compressor and pumps.

Methanol synthesis process
Frequent LC is potentially the main operation scenario involved in the flexible production of green methanol, where the flow rates of the feed gases (H 2 and CO 2 ) require adjustment. In the present study, LC operations between full-and half-load were conducted with a ramp rate (R) of 50% load per hour (considered as the condition for a basic case) and a total on-stream operation time of 15 h with load up and down.
The results of the dynamic simulations of the MS process are shown in Fig. 7. The flow rates of the feed streams (H 2 and CO 2 ) and product streams (CRD and PUR1) for the LC operation conditions are shown in Fig. 7(a). The operation load linearly decreased during t = 1-2 h and increased during t = 8-9 h. The trends of the flow rate change for CRD and PUR1 agreed with those of the LC. The main aim of PUR1 is to remove the redundant reactants and/or possible inert gases accumulated in the process (with the recycle stream); for example, the CO 2 source based on biogas production could contain inert gases such as nitrogen and methane. The gas in PUR1 contains a high concentration of H 2 (>85 mol%) and is assumed to be a fuel product in this study, which could be further utilized for other purposes such as combustion. Fig. 7(b) presents the gas composition at the reactor inlet (stream S21) for the LC operation conditions. As mentioned above, there was a high H 2 concentration in the system, and the remaining gases were primarily carbon-based (CO and CO 2 ). As shown in Fig. 7(b), after the load reduction during t = 1-2 h, slightly increasing and decreasing trends were observed during t = 2-8 h for hydrogen and carbon, respectively. Opposite trends were observed during t = 9-15 h after the load increased with faster change rates of the gas composition. The slow change rates of the gas composition in the system were attributed to the large concentration inertia of high-pressure gas present in the system (e.g., in the gas-liquid separator) compared with the small change in the gas composition. It takes more than 20 h for the gas composition to achieve a steady state (shown in Fig. S3, Supplementary data), but its influence on the system energy efficiency can be neglected within a shorter duration for steady operation (e.g., 2 h) after the LC operations (shown in Fig. 7d), and a quasi-steady-state can be assumed. Additionally, the half-load operation exhibited a slightly higher one-pass conversion of (CO + CO 2 ) in the MS reactor (R1) than that in the full-load operation; an increase from 15.7% at t = 1 h to 16.5% at t = 7 h was noted. The slight difference of one-pass conversion could be attributed to the lower pressure drop (ΔP = 1.6 and 0.4 bar at full-load and half-load, respectively) through the methanol reactor at half-load operation, which resulted in a higher average operating pressure and is favored by the methanol synthesis reactions. Fig. 7(c) shows the duty changes of the heat exchangers (HX1, HX2, and CL), and the power required of the recycle compressor (CMP) for the LC operation conditions. The duties of HX1 and CL, and the CMP power demonstrated approximately linear changes during the LC operations, which also corresponds to the trend of the feed gases shown in Fig. 7(a). The inlet cold stream and outlet hot stream of the MS reactor were well integrated by HX1 with duties of 6.21 MW at full-load and 3.10 MW at half-load. The stream S32 after HX1 (hot stream side) was cooled to 79.2 ℃ at full-load and 74.9 ℃ at half-load. However, the molar vapor fraction of S32 remains >99%; the cooler CL is necessary for further condensing the methanol and water in the gas phase of S32 with the cooling duties of 3.29 MW and 1.40 MW at full-and half-load, respectively. Additionally, the cooling of the stream S32 with low-grade heat is the main energy loss for the MS process. The objective of the heat exchanger HX2 connected to the reactor R1 was to partially exchange the high-grade heat from the hot stream S3 to the reboiler in the MD process. The heat duty of HX2 displayed small fluctuations at the minute level during the LC operations, which are attributed to the thermal inertia of the reactor (for instance, the reactor tube with stainless steel and solid catalyst) and the temperature change at the reactor (R1) outlet. Fig. 7(d) shows the energy efficiency for the MS process during the LC operations with different ramp rates (ramp time of T = 0.25-2 h). The energy efficiency for the MS process slightly increased with a decrease in load (from 87.7% at t = 0 h to 90.2% at t = 7 h), which is mainly attributed to a lower mean temperature difference (MTD) achieved for HX1 at half-load. The design of HX1 is based on the condition of 100% load with a MTD value of 24.0 ℃. The MTD value was lower (19.0 ℃) for the half-load operation, which indicated that a lower inlet temperature (hot side) for HX1 can be set to preheat the cold stream (S1) to 220 ℃. The inlet temperature for HX1 was adjusted from 240 ℃ at full-load to 232 ℃ at half-load, which resulted in a lower outlet temperature (hot side), and relatively lower energy loss through the cooling process in CL. In addition, conditions with different ramp time (T = 0.25-2 h) for the LC operations were also investigated, which displayed similar trends. Larger fluctuations were observed for the conditions with a smaller ramp time, for instance, an efficiency range of 79.3-94.1% with T = 0.25 h (corresponding to R = 200% load / h), which is attributed to the fluctuations of HX2 as mentioned above for Fig. 7(c).

Methanol distillation process
The dynamic simulation results for the MD process following MS are shown in Fig. 8. The flow rate changes of the product streams (MOH and WST) during the LC operations are shown in Fig. 8(a), and the flow rate of the feed stream CO 2 is also shown for comparison. The trends of the  product streams were close to that of the feed stream CO 2 , which changed linearly during the LC operations. Small lags were observed for the product streams after the LC operations t = 2 h and t = 9 h. The total CO 2 conversion for the MSD process was found to be 97.6% at both fullload and half-load. Fig. 8(b) presents the purity changes of the MOH and WST from the MD process and the change in the reflux ratio (flow rate ratio of REF/ MOH) for D1. The purity requirements of 99.7 wt% for the MOH was set, assuming fuel-grade methanol can be obtained for fuel blending. For example, a water content < 0.5 wt% is required in the methanol fuel standard ASTM 5797 [47]. A high mass fraction of water (>99.9 wt%) was also set, to avoid wastage of MOH from the bottom of the distillation column and to reduce the load on possible downstream water treatment. As shown in Fig. 8(b), the product purities (for MOH and WST) were well-controlled with only small fluctuations during the LC operations. The reflux ratio for the distillation column linearly increased from approximate 1.2 at full-load to 2.2 at half-load (adjusted by RAMP 3 in Fig. 2), which can be primarily attributed to the hydraulic concern of the distillation column. Compared with the condition at full-load operation, the vapor and liquid flow in the column were much lower at half-load; however, lower limits for the vapor (represented by pressure drop) and liquid flow were required to maintain efficient heat and mass transfer in the column. Consequently, a higher reflux ratio was required when the vapor or liquid flow was close to the minimum limit. In this study, commercial SULZER CY gauze packing, typically used for the MD process, was employed for the hydraulic evaluations of the distillation column (D1). The higher reflux ratio at half-load also resulted in a relatively higher reboiler duty (as shown in Fig. 8(c)), which was 1.40 MW compared with the full-load value of 1.91 MW. Additionally, under both full-load and half-load operations, 36% of the heat duty of the reboiler could be covered by the heat exchanged from HX2 (shown in Fig. 7(c)) in the MS section. Fig. 8(d) presents the energy efficiencies for the MD process during the LC operations with a ramp time of T = 0.25-2 h. Compared with the MS process, the MD process, as a traditional energy consumption process, exhibited lower energy efficiencies; moreover, the energy efficiency was lower at half-load (decreasing from 86.8% at t = 0 h to 82.4% at t = 7 h) due to the relatively higher reboiler duty required as mentioned above. Similar to the trends of the MS section, larger fluctuations of energy efficiency were observed for shorter ramp durations, for instance, the efficiency range of 78.9-92.3% with T = 0.25 h. After the LC operations, it took approximately 1-2 h for the methanol distillation system to reach a quasi-steady-state.

Methanol synthesis and distillation process
Based on the simulation results for the MS and MD processes, the total energy efficiencies of the MSD process are presented as shown in Fig. 9(b) for the LC operations with T = 0.25-2 h. The trends of the energy changes are similar to those for the MD process because of the larger fluctuations compared with those for the MS process. Storage of crude methanol (CRD) can be considered to decouple the MS and MD processes to some extent, which can improve the efficiency of the MSD process but with additional capital and operating expenses. For the basic case of T = 1 h, the efficiencies for the MSD process were 77.1% and 75.4% at full-and half-scales, respectively, corresponding to input  Fig. 9(a), and a quasisteady-state regarding energy efficiency of the whole system was achieved within 1 h after each LC operation.

Design of computer experiments
The surrogate modeling for the MSD process with LC operations was based on the dynamic model developed in Section 3.1. The KPIs for the process selected as the input (i.e., hydrogen feed flow rate, which represents the operating load) and output variables of the surrogate model are listed in Table 5. A ramp rate of R = 50% / h was set for each LC operation. The computer experiment was designed, which generated data through dynamic operations at different operating loads (full-and half-load), as shown in Fig. 10(a). The duration of steady operation after each LC operation is 2 h, assuming that a quasi-steady-state was approached based on the simulation results in Section 3.1. The total duration for the dynamic simulation, containing 34 LC operations and 32,406 data points (including the input and output variables), was 108 h with a sampling time of 72 s. The data from the operating time of t = 0-64.2 h (including 25 LC operations), t = 64.2-86 h (including 9 LC operations), and t = 86-108 h (including 9 LC operations) were used for the training, validation, and testing datasets (herein referred to as condition 1), respectively. Additionally, a testing dataset (condition 2) generated under continuous LC operations, as shown in Fig. 10(b), was also considered.

Surrogate modeling using the NARX model
The surrogate modeling was conducted for the MSD process based on the datasets generated in Section 3.2.1, and the NARX model described in Section 2.2. The values of 2 and 6 were selected for the input and feedback delays and size of the hidden neurons for the NARX model, respectively, considering the balance between model accuracy and computational time. The NARX model was trained, validated, and tested using Matlab. The modeling results exhibited good accuracy by the NARX network, and the evaluated performance for the testing data (condition 1) showed values of RMSE = 3.09 × 10 -5 , MAE = 2.30 × 10 -4 , and R 2 = 1.0. Fig. 11(a) provides a comparison of the estimated energy efficiency by the NARX model against the untrained testing dataset (condition 1), where good accuracy was achieved. Similar trends were also observed for other output variables, as shown in Table 5.
Furthermore, the constructed NARX network was tested for the extended condition with continuous LC operations, under which there was no time for the system to achieve a quasi-steady-state after each LC operation. The testing results shown in Fig. 11(b) demonstrate that the constructed NARX model is possible to track the untrained testing data   under condition 2 with relatively acceptable performance (RMSE = 0.075, MAE = 0.015, and R 2 = 0.976). Although condition 2 may not be the regular operating condition for the MSD process, there is potential for the constructed NARX network to capture its dynamic characteristics with good accuracy by adding additional training data under this condition.

Conclusions
In the present study, the MSD processes for a pilot-scale green methanol production (corresponding to 22,236 tons/year) were investigated by dynamic modeling and simulations. The LC operations for the MSD process were investigated, focusing on the energy analysis and dynamic characteristics of the system. Based on the constructed dynamic model, surrogate modeling for the MSD process was conducted using the NARX model.
For the MS process, the energy efficiency was between 87.7% and 90.2% during the LC operations (between full-and half-load) with a ramp time of T = 1 h (ramp rate of R = 50% / h). Relatively lower energy loss was achieved through the cooling process for the half-load operation, owing to the lower MTD value of the main heat exchanger HX1. Larger fluctuations were observed for the conditions with a smaller ramp time (e.g., energy efficiency with a range of 79.3-94.1% under the condition of T = 0.25 h), which was primarily attributed to the thermal inertia of the reactor.
For the MD process, the product purities for the MOH and WST can be well-controlled during the LC operations by adjusting the reflux ratio. Compared with MS, the MD process exhibited lower energy efficiencies of 86.8% at full-load and 82.4% at half-load; moreover, the energy efficiency was lower at half-load due to the relatively higher reboiler duty required; additionally, larger fluctuations of energy efficiencies (efficiency range of 78.9-92.3%) were observed for shorter ramp durations (T = 0.25 h). The total efficiencies for the MSD process under the condition of T = 1 h were 77.1% and 75.4% at the full-and half-scales, respectively.
The surrogate modeling for the MSD process by the NARX network exhibited good accuracy with the RMSE = 3.09 × 10 -5 , MAE = 2.30 × 10 -4 , and R 2 = 1.0. Additionally, the developed NARX model can track the untrained testing data for an extended condition with continuous LC operations with relatively acceptable performance (RMSE = 0.075, MAE = 0.015, and R 2 = 0.976), which can be further improved. The constructed NARX model can be further integrated with the models for other sections of the power-to-methanol process. Similarly, the applications of surrogate modeling also have potential for other power-to-x processes.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.