Bayesian inference of aqueous mineral carbonation kinetics for

Subscriber access provided by UNIV OF LOUISIANA

Process Systems Engineering

Bayesian inference of aqueous mineral carbonation kinetics for carbon capture and utilization Jonggeol Na, Seongeon Park, Ji Hyun Bak, Minjun Kim, Dongwoo Lee, Yunsung Yoo, Injun Kim, Jinwon Park, Ung Lee, and Jong Min Lee Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.9b01062 • Publication Date (Web): 15 Apr 2019 Downloaded from http://pubs.acs.org on April 17, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Bayesian inference of aqueous mineral carbonation kinetics for carbon capture and utilization Jonggeol Na,†,k Seongeon Park,‡,k Ji Hyun Bak,¶ Minjun Kim,‡ Dongwoo Lee,‡ Yunsung Yoo,§ Injun Kim,§ Jinwon Park,§ Ung Lee,† and Jong Min Lee∗,‡ †Clean Energy Research Center, Korea Institute of Science and Technology (KIST), 5 Hwarang-ro 14-gil, Seongbuk-gu, Seoul 02792, Republic of Korea. ‡School of Chemical and Biological Engineering, Seoul National University, Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea. ¶School of Computational Sciences, Korea Institute for Advanced Study (KIAS), 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea. §Department of Chemical and Biomolecular Engineering, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea. kContributed equally to this work E-mail: [email protected]

Abstract We develop a rigorous mathematical model of aqueous mineral carbonation kinetics for carbon capture and utilization (CCU) and estimate the parameter posterior distribution using Bayesian parameter estimation framework and lab-scale experiments. We conduct 16 experiments according to the orthogonal array design and additional one experiment for the model test. The model considers the gas-liquid mass transfer,

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

solid dissolution, ionic reactions, precipitations, and discrete events in the form of differential algebraic equations (DAEs). The Bayesian parameter estimation framework, which we distribute as a toolbox (https://github.com/jihyunbak/BayesChemEng), involves surrogate models, Markov chain Monte Carlo (MCMC) with tempering, global optimization, and various analysis tools. The obtained parameter distributions reflect the uncertain or multimodal natures of the parameters due to the incompleteness of the model and the experiments. They are used to earn stochastic model responses which show good fits with the experimental results. The fitting errors of all the 16 datasets and the unseen test set are measured to be comparable or lower than when deterministic optimization methods are used. The developed model is then applied to find out the operating conditions which increase the duration of high CO2 removal rate and the carbonate production rate. They have highly non-linear relationships with design variables such as the amounts of CaCO3 and NaOH, flue gas flow rate, and CO2 inlet concentration.

2


Page 2 of 62

Page 3 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Introduction

Carbon capture and utilization or carbon conversion and utilization (CCU) technology have recently drawn attentions in both industry and academia because of their economic and environmental potentials. 1,2 Current development progress of CO2 utilization is at the technology readiness level 6 (TRL6), which is the stage of pilot plant demonstration. 1 There are several categories of technology under CCU that can convert CO2 to useful chemicals and fuels. Some of them use electro-catalytic, photo-electrochemical, photo-catalytic, and thermochemical reduction, 3 and others use CO2 as a feedstock to produce various chemical products such as urea, methanol, and salicylic acid. 2 Mineral carbonation, one of promising CCU technologies, captures CO2 in the form of mineral carbonates that can be fabricated as building compartments after post-processing. 4,5 It is proved to be at least carbon-neutral, and carbon-negative when atmospheric CO2 is used as feedstock. The advantages of this technology are manifold. First, it can directly use the flue gas from power plants or incineration plants without any preprocessing steps. Second, it can use the industrial wastes such as steel slags, fly ash, wastewater and brine as mineral sources. 4,6 In addition, the solid carbonate products are thermodynamically stable and have a mature market as construction materials to anticipate the economic potentials. 6 Several companies have been actively engaged in the development of mineral carbonation process and made visible progress. Regarding industrial applications, an Australian based start-up, Mineral Carbonation International (MCi), invested by Orica since 2013, has operated CO2 mineral carbonation pilot plant since 2016. 7 They announced a plan to construct a full-scale production plant (20,000-50,000 tons/year of carbonate and silica by-products) by 2020. 8 The pilot plants of Calera Corporation capture up to 2 tonnes/day of CO2 into calcium carbonate product directly from the raw flue gas from power plants. They also constructed a commercial fiber cement board line integrated with the mineral carbonation pilot plant to produce full-size fiber cement board sheets. 9 The Carbon Capture Machine developed by the University of Aberdeen has reached TRL4, currently producing 200 kg/day 3



of precipitated calcium carbonates. COSIA Carbon XPRIZE Challenge, which was launched in 2015, is now supporting its further development to TRL6 demonstration plant. 10 Daewoo E&C in South Korea developed and operated a pilot plant which can capture 10 tonnes/day of CO2 in 2012. They built the scaled-up plant which can process 40 tonnes/day of CO2 in 2017, 11,12 in which the authors have participated. U.K. start-up Carbon8 has a facility which treats 2,000 tonnes/year of CO2 and produces building aggregates using wastes such as ashes from municipal incinerators and energy plants. 13 Although a lot of effort is being made to commercialize mineral carbonation processes as the main technology of CCU, there still remain challenges. First of all, the reaction rates are slow due to the mass transfer resistance between different phases - the system involves the dissolutions of mineral and hydroxide ions from the solid reactants, e.g. Ca(OH)2, and carbonate ions from the flue gas. In order to overcome the slow reaction rates, it is essential to optimize the reactants and the design of the reactor by modeling the entire reaction kinetics including the interphase mass transfer. The good reactor design can help accelerate the mass transfer rates by homogeneously distributing the gas bubbles and solid particles. Many studies on mineral carbonation kinetics and the reaction parameters have been reported and verified with experimental results. 14–19 All of them used deterministic methods for the parameter estimation and suggested deterministic values which best fit their experimental results. However, some of the parameters, especially the ones related to hydrodynamics and mass transfer, are likely to have different values under different operating conditions and disturbances, such as the fluctuation of the CO2 concentration in the flue gas, the inhomogeneous size distributions of solid particles. Therefore, the use of the deterministic parameters in the model may hinder the accurate predictions of output responses. Bayesian parameter estimation can find the posterior probabilistic distributions of the parameters by taking account of their distributed natures as well as the measurement errors of the experiments. The model can then be updated to a stochastic model using uncertainty quantification methods 20,21 and used for stochastic model predictive control (SMPC). 22,23 4


Page 4 of 62

Page 5 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Although there is no general Bayesian parameter estimation framework available, we can find a lot of applications in many different fields including chemical engineering, 24,25 carbon mitigation, 26,27 bioengineering, 28 cognitive science, 29 cosmology, 30 electrical/electronics engineering, 31 and artificial intelligence. 32 One challenge of Bayesian estimation is that a huge amount of samples are required. Thus, when the computational load of the model is too expensive, function approximation methods such as second order polynomial 25,33 and polynomial chaos expansion (PCE) 24,34 are adopted together. Another big challenge is dealing with the multiple datasets with different levels of uncertainties. Experiments are usually performed under varing conditions (or design variables), so the amount of data and uncetainty inevitably varies as well. Therefore, a proper Bayesian parameter estimation framework which can handle these problems is demanded to resolve the issues in mineral carbonation technology. Here, we propose fundamental aqueous mineral carbonation kinetics model for threephase (solid: Ca(OH)2, NaOH, and CaCO3; liquid: absorbent solution with ions; gas: flue gas) reactor. The model takes a form of differential algebraic equations (DAEs). Using an existing DAE solver, 35 we solve stiff and hybrid systems with significant magnitude differences among the parameters and discrete events such as the disappearance of solids after complete dissolutions. The parameters in the model are estimated through a Bayesian parameter estimation framework, which can quantify uncertainties in the experimental outputs. In order to show the efficacy of our framework and obtain the posterior distributions of mineral carbonation kinetic parameters, we constructed a lab-scale process and conducted the experiments according to an orthogonal array design of experiments (DoE). This paper is organized as follows. In Section 2, we present the experimental methods of lab-scale mineral carbonation process. In Section 3, we describe the mathematical model of reaction kinetics for 3-phase mineral carbonation reactor. A Bayesian parameter estimation scheme is explained in the next Section 4. In Section 5, we discuss physical interpretation of the experimental results and the parameter estimation. Finally, Section 6 provides 5



concluding remarks.

6


Page 6 of 62

Page 7 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2

Experimental Methods

2.1

Solution and gas preparation

The absorbent solution was a mixture of Ca(OH)2 (purity > 95.0% Ca(OH)2 ACS reaction, Sigma Aldrich) and NaOH (purity > 97.0% NaOH ACS reaction, Sigma Aldrich) in distilled water. CO2 and N2 gas mixture was used as an artificial flue gas. The CO2 volume ˙ concentration of inlet flue gas, φinlet CO2 , and the flue gas flow rate, V were controlled by mass flow controller (MFC) at 101.325 kPa. The reactor temperature, T , the weight fractions of ˙ Ca(OH)2 (wCa(OH)2 ) and NaOH (wNaOH ), φinlet CO2 , and V were decided by an orthogonal array experiment design method (a strength 2, 4 levels, 5 factors and index 1; 2-(4,5,1)). We conducted 16 experiments according to the experiment design and one additional experiment for validation, as described in Table 1.

2.2

Laboratory-scale mineral carbonation process

The overall schematic process flow diagram is shown in Figure 1. 500 ml PYREX® doublejacketed reactor was used as a gas mixer and a carbonation reactor. To filtrate the solid particles including calcium carbonate, calcium hydroxide, and sodium hydroxide, glass fiber filter with 1.6 µm pore was applied at the inlet blower of CO2 gas. The temperature of the gas mixer and the carbonation reactor, T , were controlled by a thermostat water bath. To prevent concentration changes caused by evaporation of the absorbent solution, a condenser was installed with a water chiller. Experiment outputs were CO2 volume concentration of the outlet flue gas, φoutlet CO2 , and pH of the absorbent solution in the reactor. They were recorded from the gas analyzer (multi-gas analyzer, Sensoronic) and the pH meter (HM-41X, DKK-TOA) in time-series (sampling time was 30 seconds), respectively. To start-up the process, temperatures of water baths of the gas mixer and the reactor were first set following the experiment design. In order to get perfectly mixed flue gas as well as to prevent fluctuations of the concentration and the temperature, the CO2 gas was sent to a 7



vent line for a while until the temperature reached a steady-state. Meanwhile, highly pure N2 gas was purged to the reactor to remove the latent dissolved air and impurities. Afterwards, the N2 purging was stopped and the reactor was filled with 500 mL of the absorbent solution. Lastly, we opened the valve prior to the reactor and made the flue gas blown into the reactor. The reacting mixture was stirred in 300 rpm using a magnetic bar to distribute bubbles and solid particles. The experimental outputs were recorded every half a minute. When the CO2 volume concentration of the outlet flue gas showed no difference to the inlet gas, which means that the absorbent solution cannot absorb CO2 gas any longer, we assumed that the absorbents in the solution were entirely consumed. So we stopped logging the data and shut down the process. Similar experimental settings and procedure were applied in our previous papers as well. 36,37

8


Page 8 of 62

Page 9 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3

Mathematical Models

3.1

Kinetics of aqueous mineral carbonation process

The reactions involved in the aqueous mineral carbonation process can be categorized into four groups: First, mass transfer of carbon dioxide gas into the water solution, which is the most important reaction that decides the overall CO2 removal efficiency. Second, dissolutions of two solid species, NaOH and Ca(OH)2. Third, ionic reactions between the dissolved ions derived from the gas and solid reactants. Last, the precipitation of CaCO3, a final product. Generally, the ionic and precipitation reactions are fast, whereas the gas-liquid mass transfer and solid dissolutions are slow and regarded as rate-limiting steps because of the resistance between two different phases. 3.1.1

Mass transfer between gas and liquid phases CO2(g) → CO2(aq).

(1)

This is a physical absorption of CO2 gas to the liquid phase, whose transfer rate per unit gas volume is expressed as rmass transfer = (θkl kl )Ag (θE E)(H CO2 [CO2(g)] − [CO2(aq)]),

(2)

where kl and E are the overall mass transfer coefficient and the enhancement factor by chemical reactions, to which the adjustment factors θkl and θE are multiplied respectively, in order to compensate the uncertainties in the nominal parameters. Ag is the surface area per unit gas volume, H CO2 is the Henry constant of CO2 in electrolyte solutions, and [ · ] denotes the molar concentrations of corresponding species. The overall mass transfer coefficient, kl , is derived from Sherwood (Sh) relation, 38

Sh =

kl θdb = 2 + 0.0015Re0.89 Sc0.7 , DCO2 9


(3)


Page 10 of 62

where θdb is the bubble diameter, which we assumed to have a constant value that must be found out by the parameter estimation. DCO2 is the diffusivity of CO2 under alkaline environment, and Re and Sc are Reynolds and Schmidt numbers, respectively. The correlation of DCO2 is given by DCO2 CO Dw 2

DwCO2

= 1 − 1.29 × 10−4 [OH–], −2119 −6 , = 2.35 × 10 exp T

(4) (5)

where DwCO2 is the diffusivity of CO2 in pure water. 39 Although we can calculate the mass transfer coefficient at each time step using the above correlations, uncertainties still exist in the hydrodynamics related terms like Re and Sc. Thus, we introduce an adjustment factor, θkl , to compensate these uncertainties, and estimate it with the experimental results. The surface area per unit gas volume, Ag , is determined from the given bubble diameter, θdb , and the gas holdup, αg , as Ag =

αg . θdb

(6)

θdb and αg are uncertain but have significant influences on the mass transfer rate. Thus, it is advantageous to estimate them as parameters, rather than giving arbitrary values. In order to reflect the change of gas holdup according to the gas inlet, we introduce a correlation with three parameters, θa , θb , and θc : αg = θa + θb Ugθc ,

(7)

where Ug is the superficial gas velocity. To sum up, four parameters, θdb , θa , θb , and θc , are to be estimated.

10


Page 11 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The enhancement factor, E, is calculated from the correlation 40 of

E=

   −

Ha2 2(E∞ −1)

+

q

Ha4 4(E∞ −1)2

+

E∞ Ha2 E∞ −1

+ 1, E∞ > 1,

  1,

(8)

E∞ ≤ 1,

where r CO [OH–]DCO2 D 2 E∞ = 1+ –, 2DCO2 H CO2 [CO2(g)] DOH p k11 DCO2 [OH–] . Ha = kl

(9) (10)

–

DOH is the diffusivity of OH− ion, which takes the value of 5.3×10−9 , 41 and k11 is the rate constant of a reaction between CO2(aq) and OH− (aq), which is discussed in the next part. The adjustment factor, θE , is multiplied to the Equation (8) when calculating the mass transfer rate. It is to represent all uncertain parameters involved in the Equations (8)-(10). The Henry constant for electrolyte solution, H CO2 , is obtained from log

HwCO2 H CO2

=

X

(11)

(hi + hg )ci ,

i

HwCO2

−7

= 3.54 × 10 RT exp

2044 T

,

(12)

where HwCO2 , hi , hg and ci are the Henry constant of CO2 for pure water, the ion-specific parameter, the gas-specific parameter, and the concentration of each ion i, respectively. Equation (11), (12) and all the parameter values are adopted from literature. 42,43

11



3.1.2

Page 12 of 62

Solid dissolutions

kf

2+ – * Ca(OH)2(s) − ) − Ca (aq) + 2OH (aq), kb

kf 2

2+ – −− * CaOH+(aq) ) − − Ca (aq) + OH (aq), kb2

kf 3

+ – −− * NaOH(s) ) − − Na (aq) + OH (aq). kb3

(13) (14) (15)

The dissolutions of alkaline solids provide mineral and hydroxide ions. All the forward and backward reaction rates of these reactions can be expressed as products of each reactant’s concentration and the corresponding rate constants. The equilibrium and rate constants are given in Table S1 in Supplementary Material. It should be noted that kf and kb in Equation (13) are dependent on total surface area of the solid to reflect the decrease of dissolution rate owing to the shrinkage of the solid particles as the reaction proceeds. The value of kf 2 was randomly set to 1,000 because we could not find a relevant source. However, we did not perform parameter estimation on it because the results were not affected by its value at all. kf 3 and kb3 were also randomly set to 109 and 0, respectively. We checked that they have negligible impacts on the results as long as NaOH completely dissolves in a short time. 3.1.3

Ionic reactions

k

11 – − * CO2(aq) + OH–(aq) − ) − − HCO3 (aq),

k12

k

21 2– −− * HCO3–(aq) + OH– ) − − CO3 (aq) + H2O(l),

k22

k

31 − * H+(aq) + OH–(aq) − ) − − H2O(l),

k32

k

41 – + −− * CO2(aq) + H2O ) − − HCO3 (aq) + H (aq).

k42

(16) (17) (18) (19)

The above reactions represent the conversion of CO2(aq) into bicarbonate and carbonate ions in an acid-base equilibrium. The forward and backward reaction rates follow the elementary 12


Page 13 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reaction rate equations, like in the solid dissolutions. The equilibrium constants and rate constants are given in Table S2 where the equilibrium constants for Equations (16) -(19) are denoted as K1 , K2 , Kw , and K4 , respectively. 3.1.4

Precipitations

k

51 − * Ca2+(aq) + CO32–(aq) − ) − − CaCO3(s),

k52

k

61 + −− * Ca2+(aq) + HCO3–(aq) ) − − CaCO3(s) + H (aq).

k62

(20) (21)

Two reactions here explain production of the main product, calcium carbonate. Neither of these becomes a rate-limiting step because they are fast ionic reactions producing the thermodynamically stable product. The equilibrium and rate constants are given in Table S3.

3.2

Differential algebraic equation (DAE) model for the reactor

We constructed a set of differential algebraic equations (DAEs) based on the kinetics described in the previous parts. The exact formulations as well as the assumptions we made are explained in this section. Some algebraic equations like rate constants or equilibrium constants which are omitted here can be found in Supplementary Material. 3.2.1

Assumptions

• The reactor consists of 10 vertical compartments with a same volume where the gas passes through each of them from the bottom. The concentration of CO2(g) declines as it moves to the upper compartment (Figure 2). • The gas distribution is homogeneous everywhere, i.e. gas holdup has a constant value in all the compartments.

13



• Gas bubbles have a constant diameter and do not shrink during the reactions. It is reasonable because the portion of CO2 is small compared to N2 . • CaOH2 particles shrink as the reaction proceeds. We fix the number of particles and make the total volume, surface area, and diameter changed. When CaOH2 is completely consumed, the surface area becomes zero and the dissolution no longer takes place. • Solids and ions are uniformly distributed throughout the reactor, i.e. their concentrations have the same value in every compartment.

14


Page 14 of 62

Page 15 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2.2

Model formulations

[CO2(g)i−1 ] [CO2(g)i ] d[CO2(g)i ] = Qi−1 − Qi dt Vgas Vgas 6(θa + θb Ugθc ) −θkl kl θE E(H CO2 [CO2(g)i ] − [CO2(aq)]), i = 1, ..., N, (22) θdb d[CO2(aq)] = −k11 [CO2(aq)][OH–] + k12 [HCO3–] − k41 [CO2(aq)] + k42 [HCO3–][H+] dt N 6(θa + θb Ugθc ) Vgas X + θE E(H CO2 [CO2(g)i ] − [CO2(aq)])), (23) (θkl kl VR i=1 θdb d[HCO3–] = k11 [CO2(aq)][OH–] − k12 [HCO3–] − k21 [HCO3–][OH–] + k22 [CO32–] dt +k41 [CO2(aq)] − k42 [HCO3–][H+] − k61 [Ca2+][HCO3–] + k62 [H+], d[CO32–] = k21 [HCO3–][OH–] − k22 [CO32–] − k51 [Ca2+][CO32–] + k52 , dt d[H+] = −k31 [OH–][H+] + k32 − k42 [HCO3–][H+] + k41 [CO2(aq)] dt +k61 [Ca2+][HCO3–] − k62 [H+],

(24) (25)

(26)

d[OH–] = −k11 [CO2(aq)][OH–] + k12 [HCO3–] − k21 [HCO3–][OH–] dt +k22 [CO32–] − k31 [OH–][H+] + 2θAs As (kf − kb [Ca2+]f 4 [OH–]2 f 2 ) +kf 2 [CaOH+] − kb2 [Ca2+][OH–] + kf 3 − kb3 [Na+][OH–],

(27)

2+

d[Ca ] = −k51 [Ca2+][CO32–] + k52 − k61 [Ca2+][HCO3–] + k62 [H+] dt +θAs As (kf [Ca(OH)2] − kb [Ca2+]f 2 [OH–]2 f 4 ) +kf 2 [CaOH+] − kb2 [Ca2+][OH–] + kf 3 ,

(28)

+

d[CaOH ] = −kf 2 [CaOH+] + kb2 [Ca2+][OH–] + kf 3 , dt d[Na+] = kf 3 − kb3 [Na+][OH–], dt d[CaCO3] = k51 [Ca2+][CO32–] − k52 + k61 [Ca2+][HCO3–] − k62 [H+], dt

15


(29) (30) (31)


Qi = Qi−1 −

6(θa + θb Ugθc ) RT θkl kl θE E(H CO2 [CO2(g)i ] P θdb (32)

−[CO2(aq)])Vgas , VCaOH2 = As = ns = log(f ) = I =

1

Page 16 of 62

0 (ZCaOH − ([Ca2+] + [CaOH+] + [CaCO3])M WCaOH2 VR ), 2

ρCaOH2 2 π 1 ( ) 3 (6VCaOH2 ) 3 , ns Z0 , θd0 4 s 3 ρ π( ) CaOH 3 √2 2 −0.5 I √ , 1 + 1.4 I 0.5(4[Ca2+] + [H+] + [OH–] + [HCO3–] +4[CO32–] + [CaOH+] + [Na+]).

(33) (34) (35) (36)

(37)

In Equation (22), [CO2(g)i ], Vgas , and Qi are the concentrations of CO2(g)i in the ith compartment, the volume of mixed gas in one compartment, and the volumetric gas flow rate from the compartment i-1 to i, respectively. [CO2(g)0 ] is the inlet CO2(g) concentration, 0 which can be calculated from the design variable, φinlet CO2 , using the ideal gas law. Q is the

volumetric gas flow rate, which is equivalent to the design variable, V . Vgas can be calculated as the product of the gas holdup and the compartment volume, i.e., Vgas = (θa + θb Ugθc ) × VNR . VR in Equation (23) is the volume of the reactor. As in Equation (27), is the total surface area of CaOH2 particles, θAs is an adjustment factor to compensate the uncertainty in As , 0 and f is the activity coefficient. In Equation (33), ρCaOH2 , ZCaOH , and M WCaOH2 are the 2

density, molecular weight and initial weight of CaOH2 , respectively. In Equations (34) and (35), ns is the number of the particles and θd0s is the initial diameter of CaOH2 particles, which is one of the important but uncertain parameters that we perform parameter estimation. The other algebraic equations such as k11 , H CO2 , kl , and E are already explained in previous sections.

16


Page 17 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2.3

Parameters (θ) to estimate

To sum up, there are eight parameters (θ) to estimate in total: θa , θb , and θc are the parameters which are related to the gas holdup, θdb and θd0s are the initial diameters of gas bubbles and solid particles, and θkl , θE , and θAs are the adjustment factors for uncertain parameters in the reactor model. 3.2.4

Output responses (y)

Two output responses, y1 (CO2 volume concentration of outlet flue gas) and y2 (pH), were monitored during the experiments. From the proposed model, they are defined as

3.3

y1 = φoutlet CO2 ,

(38)

y2 = − log(Kw ) + log([OH–]).

(39)

Discrete events for simulation procedure

Three discrete events, that we need to reinitialize the DAE model, can happen. The first is the entrance of the gas to the reactor. In the experiments, the gas is allowed to enter the reactor after given a time for NaOH and CaOH2 to sufficiently dissolve in water. We reproduced the same procedure in our model by first running the DAE model without CO2related variables. The second and third are the complete dissolution of NaOH and Ca(OH)2, respectively. After the complete dissolution of solid particles, the dissolution kinetics of each solid species were deactivated. If not, Na+, Ca2+, and OH– ion can be generated infinitely even after the sources are gone. Furthermore, if we just force kinetic parameters to be zero rather than finding an exact event time using a root finding algorithm, convergence problems can happen. To prevent this, we introduced an external functions, g([Ca(OH)2])

17



Page 18 of 62

and g([NaOH]), that calculate the remaining weights of Ca(OH)2 and NaOH: 0 g([Ca(OH)2]) = ZCaOH − ([Ca2+] + [CaCO3] + [CaOH+])M WCa(OH)2 , 0 g([NaOH2 ]) = ZNaOH − [Na+]M WNaOH ,

(40) (41)

where Z 0 and M W are the initial weight and molecular weight of each species, respectively. When these functions give zero values, kf , kb , kf 3 , and kb3 are set to zero and the model is reinitialized.

3.4

Numerical setting

The proposed mathematical model was written in Matlab script. We used the Sundials IDAS solver with the dense linear solver and zero-crossing root-finding algorithm 35 to handle the stiff and hybrid DAEs system. Furthermore, since the kinetics show stiff behaviors like 108 109 order-of-magnitude difference among parameters, the model occasionally fails to converge within a feasible time. Thus, we monitored the divergence of the solver and tuned solver parameters such as relative tolerance (10−7 ) and absolute tolerance (10−7 ) to make as many samples converged as possible. We simulated our model using Intel Xeon E5-2667 v4 (3.2 GHz) and the average simulation time for one scenario was less than one second.

18


Page 19 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4

Bayesian Parameter Estimation

One objective of our study is to infer posterior distributions P(θ) of the unknown parameters, θ, from experimental data. Our problem has characteristics that data acquisition is expensive, a non-adpative black-box (a factorial method) 44–47 is used as an experiment design method, and the mathematical model is complex to perform Markov chain Monte Carlo (MCMC) sampling. Taken together, we are tackling the Bayesian parameter estimation problem with a small number of data sets and complex model.

4.1

Problem formulation

Let x be a set of design variables for an experiment, and y the observed response. Specifically, we consider the case where the response is a time series: y = {y1 , y2 , · · · , yt , · · · , yT }, written in the form of a T -dimensional vector. Each experiment reports multiple response types, {y1 , · · · , yD }, where each response yd is a T -dimensional vector. We have D = 2 in our problem: y1 measures the CO2 volume concentration of outlet flue gas, and y2 the pH of the solution. In our experiments, the response length T is fixed for y1 and y2 from the same experiment, because they are measured simultaneously (Figure 3). We consider a finite sets of design variables {x(1) , · · · , x(M ) } that correspond to multiple experiments, where ˙ x = (T, φinlet CO2 , wCa(OH)2 , wNaOH , V ). We have M = 16 experiments; M cannot be very large in most cases, because experiments are costly. The length of response T may be different for different experiments. With M = 16 experiments and D = 2 response types, we have (m)

M × D = 32 datasets; by a single “dataset” we mean {x(m) , yd } of a particular (m, d), where m = 1, · · · , M and d = 1, · · · , D. Let us denote the first-principle model in previous section by f . Given a set of design variables x and a fixed set of parameters θ, the model predicts a response time series as ymodel = f (x, θ). Here θ is a K-dimensional parameter vectors, where K = 8 in the current problem. Specifically, we have θ = (θE , θAs , θa , θb , θc , θdb , θd0s , θkl ). In general, there is always

19



Page 20 of 62

some discrepancy (or model-plant mismatch) between the actually observed response y and the modeled response ymodel . We can write this as y = f (x, θ)+, where is the error vector. We assume that the error is distributed according to a multivariate normal distribution, with a T × T covariance matrix Σ: 25,33

= y − f (x, θ) ∼ N (0, Σ).

(42)

For a complicated kinetic model, calculation of the model response f (x, θ) while varying the value of θ can be computationally demanding. To make computations affordable, various surrogate models can be applied such as neural networks, Gaussian process, and PCE. In this study, we approximate the modeled response using a quadratic hyper-surface to avoid additional training step: f (x, θ) ≈ ˜f (θ) = c + b> θ + θ > Aθ,

(43)

where ˜f (θ) is the surrogate model response, A is a symmetric K × K matrix, b is a vector of length K, and c is a scalar, and K is the dimensionality of θ. The parameters were sampled on a Central Composite Design (CCD) and a Latin hypercube sampling (LHS) 48 (1,000 samples in total) within the prior hypercube. The quadratic surface fit was performed separately for each dataset (m, d), to obtain an approximate response function ˜f (m,d) such (m)

that yd

4.2

≈ ˜f (m,d) (x(m) , θ).

Bayesian posterior inference

According to the Bayes rule, the posterior distribution P(θ) is formulated like

P(θ) ≡ p(θ|data) ∝ p(data|θ) · p(θ),

20


(44)

Page 21 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where p(data|θ) is the likelihood of observing the data from a model parameterized by θ, and p(θ) is the prior distribution on θ. In the rest of this section, we outline how we can construct the posterior distribution by calculating the likelihood and specifying a prior distribution, and how the posterior can be sampled by the MCMC method. 4.2.1

Likelihood

According to the Gaussian error assumption, the likelihood of a given set of data depends not only on θ, but also on the covariance matrix Σ. In this sense, Σ is another parameter in the inference problem, although it is not of our primary interest; we are only interested in the distribution of θ. A full-fledged Bayesian approach for this problem is to consider the joint parameter space of (θ, Σ) and then to marginalize over Σ at the end, as in some of the previous works. 25,33 Here we take a different approach, in which we estimate the error ˆ directly from the data and an inferred distribution over θ, and simply assume covariance Σ ˆ In other words, that the distribution over Σ is highly localized around its “true” value, Σ. we aim to calculate Z p(data|θ) =

ˆ dΣ p(data|θ, Σ) p(Σ) ≈ p(data|θ, Σ).

(45)

ˆ is estimated separately for each dataset, reflecting the fact that each The error covariance Σ observation may come with a different amount of uncertainty. Specifically, we parameterize the covariance matrix in an exponential form, so that the covariance element between the two time-points ti and tj is modeled as Σij = σ 2 exp(− |ti − tj | /τ ), and estimate the fluctuation scale σ and the de-correlation timescale τ . Estimation of the fluctuation scale σ is performed by taking an average over the posterior distribution, which requires an iterative algorithm. See Figure 4 for a schematic illustration; also see Supplementary Material for more details. For multiple datasets, in general, the likelihood p(data|θ) involves a multiple integral over

21



Page 22 of 62

the Σ’s; if we consider all M experiments and D responses, the total likelihood is written as Z p(data|θ) =

Z dΣ1,1 · · ·

dΣM,D p(data|θ, Σ1,1 , · · · , ΣM,D ) p(Σ1,1 ) · · · p(ΣM,D ).

(46)

With the localization assumption of Equation (45), the total likelihood can be simply separated as p(data|θ) ≈

M Y D Y

(m)

p((x(m) , yd )|θ, Σm,d ).

(47)

m=1 d=1

Therefore, the joint likelihood of multiple datasets can be expressed as the product of the individual-dataset likelihoods. 4.2.2

Prior

We used a uniform prior distribution with a range constraint, such that p(θ) ∝ 1 if θ ∈ C and p(θ) = 0 otherwise, defined by a hypercube C ⊂ R8 . The detailed boundary values of the hypercube C , which are basically ranges of each parameter, are described in Table 2. In order to ensure that the entire parameter space is efficiently explored by the MCMC sampling, we considered the log of a parameter elements when the prior range spans multiple orders of magnitude.

4.3

Sampling

We performed MCMC sampling with the Metropolis-Hastings algorithm 49 to approximate the posterior distribution. We used an iterative sampling algorithm to adaptively adjust the proposal distribution of the Metropolis-Hastings sampler, as well as to estimate the error covariance at the same time. See Supplementary Material for the details. Once we have sampled a chain of parameters, {θ 1 , · · · , θ N }, the response to a given design variable x that is predicted by the posterior is:

fmean (x) =

N 1 X f (x, θ n ), N n=1

22


(48)

Page 23 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where N is the number of parameter samples in the chain. Similarly, we can consider the mean of all obtained samples as the representative value of the posterior distribution:

θ mean

N 1 X θn . = N n=1

(49)

It is also possible to extract mode (θ mode ) among a chain of parameters in two-dimensional parameter space. However, it is difficult to extract mode from the multi-dimensional parameter space because of the difficulty of defining the unit volume for calculating the number samples. Thus, we sequentially extracted the mode from two-dimensional parameter space and used the final one.

23



5

Results and Discussion

5.1

Stochastic output response

When the entire process of Bayesian parameter estimation is completed, the probability distributions of the output responses as well as θ can be obtained as in Figure 5. It shows that the deterministic model can be non-intrusively extended to the stochastic model through Bayesian parameter estimation. The uncertainties from the model-plant mismatch or the plant disturbance are reflected in the parameters and propagated to the output responses through them. The propagated uncertainties in the output responses are visualized in Figure 5a with the standard deviation or the confidence interval (CI). However, there are occasions where we need to pick a single response curve rather than handling a stochastic model. A straightforward way is to choose a single representative value of θ like θ mode or θ mean from the posterior distribution and solve the model with it to earn a single response, i.e. f (θ mode ) or f (θ mean ). Another way is to sample many θ from the posterior, solve the model using each of them and take the average to get the representative response, i.e. fmean . If the posterior had the normal distribution, all three of them, namely f (θ mode ), f (θ mean ) and fmean , would follow an exactly same trajectory. However, when the posterior has multiple peaks like in Figure 6, it is beneficial to adopt fmean among them because f (θ mode ) has a risk that θ mode is not picked from the highest peak and f (θ mean ) could be inferior to f (θ mode ) when f is not monotonic with respect to θ. In Figure 5a, where we show f (θ mode ), f (θ mean ) and fmean , it is noticeable that f (θ mean ) and fmean are almost overlapped, suggesting that the response is changing monotonically with respect to θ. Based on these observations, we concluded to use f (θ mean ) and fmean throughout the further discussions when we need a deterministic single response curve.

24


Page 24 of 62

Page 25 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5.2

Quality of parameter estimates

We were able to fit most of the observed responses with good accuracies (Figure S1-S3). This means that our DAE model captures the important features of a complicated CCU system in which multiple phases and reactions are involved in a discrete manner and that our parameter estimation method works. The goodness of fit of the model responses is quantitatively measured and compared with the results from several optimization methods using the following fitting error equation, v u M D uX X (y(m) − ˜f (m,d) (x(m) , θ))(y(m) − ˜f (m,d) (x(m) , θ))> d d . error = t 2 (m) (m,d) (σ ) T m=1 d=1

(50)

We used three different solvers, a genetic algorithm (GA), a dividing hyper-rectangle (DIRECT), 50–52 and a sequential quadratic programming (SQP), to find the θ opt which minimizes Equation (50). We set the maximum number of iterations as same as in the Bayesian method for a fair comparison. In Figure 7, the fitting errors of f (θ mean ), f (θ mode ), fmean and f (θ opt ) for all 17 experiments and 2 responses are shown. Overall, none of them shows dominant performance to others. Although the variance of errors is higher in the Bayesian estimation than in the optimization, the average of them have similar magnitudes in both methods. It means that our method can provide deterministic parameter values which are as feasible as optimizations can give while suggesting the probability distributions as well. In addition, we left one dataset unseen during the parameter estimation and measured the fitting errors when the obtained parameters were applied. The resulting test set errors are acceptable compared to the other experiments. Particularly, the errors from our method are distinctively lower than the errors from the optimizations in both responses. It supposes that our method is better in covering the unseen operating ranges than deterministic approaches are. However, Exp. 1 and 14 seem to show particularly higher errors than the others. We suspect that it is mainly due to two reasons. First, since the MCMC sampling cannot perform 25



sampling outside the parameter boundaries, the optimal parameter cannot be estimated corinlet rectly near the boundary. Second, the absence of NaOH makes φoutlet CO2 reach φCO2 very quickly

and the experiment terminate early. Thus, T(m) is shortened and the error is increased. We also devide the error by T(m) in Equation (50). In addition, short experiments cannot reflect the slow dynamics of CaOH2. However, Exp. 7 and 12 show lower errors than Exp. 14 although they are lack of NaOH as well. In Exp. 14, the removal efficiency converges to 0 and the experiment is over before pH drops enough. We suppose that our algorithm, which focuses on stiff response dynamics through de-correlation timescale τ , recognized little contribution on this dataset compared to the others.

5.3

Assessment of parameter uncertainties

Final MCMC sampling was conducted with 5,000 samples, 10 maximum iterations, tight chain criterion, and 0.5 step scale factor at each iteration step. At the final sampling, the number of samples and burn set to 20,000 and 2,000, respectively. More details about the numerical settings are provided in Supplementary Material. In Figure 8, we visualized two-dimensional joint marginal posterior distributions as well as marginal posterior distributions of the parameters which were sampled using the entire datasets (all 16 experiments and 2 responses). In Figure 9, we provide each parameter’s marginal posterior distributions inferred from each experiment. We can clearly see that every experiment results in a different posterior distribution. This is due to the errors derived from both model and experiments. If the model were perfect and the parameters perfectly independent of the design parameters, all the experiments would bring a same posterior distribution. However, there are some dependencies between the parameters and the design variables, that we could not include in the model like the effect of suspending solids on the gas holdup and the reaction rates.It means that even the true values of the parameters can be varied under different conditions. In addition, the potential amount of error of each experiment would be all different. Some experiments were more sensitive to 26


Page 26 of 62

Page 27 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


external disturbances than the others. For example, an experiment with higher solid loading has more chance to have the gas sparger clogged, or the one with lower gas volumetric flow is more susceptible to measurement errors. As we investigate characteristics of the posterior distributions in Figure 8 and Figure 9, we are able to categorize the parameters into three groups. The first group consists of θE and θkl , whose posteriors show unimodal behavior. In this group, the posteriors inferred from all experiments are narrowed to the area where the posterior distributions from each experiment are mostly concentrated. These parameters are likely to have deterministic true values derived from physical and chemical natures of the system. For example, θE has a value greater than 1 because the pH could have a greater impact on the response in reality than in the model, and θkl has a value less than 1 probably because of the existence of the suspensions which hinder the mass transfer between gas and liquid. The second group are θa , θb , θc and θdb , which have broad distributions not only in the single-experiment posteriors but also in all experiments posterior. It is interesting that the bubble diameter, θdb , turns out to show a distribution although we regarded it as a constant when we first designed the model. The third group includes θAs and θds , which show multimodal behavior in their posteriors. We suppose these two parameters might have inherent multimodal natures because their posteriors from single experiments are pointing at very different values. However, the posterior from all experiments shows only one sharp peak. This is because more local peaks are observed when more data are used, but MCMC algorithm cannot easily escape from a local peak once it is trapped. We can partly explain the multimodal or broad natures of the second and third parameter groups with their many-to-one properties; we can find multiple sets of (θAs , θds ) and (θa , θb , θc ) that produce very similar results. Also, the core-shell structure of the solid particle, which is formed by the precipitation of CaCO3 on the surface of CaOH2, 53 can be the reason as well because it must show completely different dissolution kinetics than original CaOH2. However, multimodal distributions are hardly observed when multiple datasets are used 27



for the inference. In this case, we can tell the existence of other local peaks by increasing the tempering factor. Sampling with a larger tempering factor means accepting worse solutions with a higher probability; graphically, it is equivalent to sampling from a flatter (or tempered) distribution than the original one which has a tempering factor of 1. Although the tempered distribution is not as same as the original distribution, it is useful when learning the overall shape of the distribution, especially when detecting small local peaks. However, the tempering factor does not affect the position of the biggest peak which has a dominant effect on the mean and mode values of the parameters. Figure 9 (see dashed lines) and Figure S3 clearly show that the increase of the tempering factor can reveal broad or multimodal natures of the posterior distributions. For more investigations, we track how θ mode , θ mean and the joint posterior are changed as the tempering factor increases (Figure 6 and Figure 10). θ mode shows a significant movement as no distinctive global peak exists. θ mean is close to θ mode at low tempering factors where only one peak is detected, but moves far from θ mode to middle of newly observed peaks when the tempering factor increases. θ mode would give a better fit for a specific response although it could be too sensitive sometimes, whereas θ mean would give more reasonable response when there are multiple datasets and the best parameters found out from each dataset is different. We found out that a strong correlation exists between θAs and θds from Figure 10. We suppose that the model underestimates the solid dissolution rate when the particles are large and the other way around when the particles are small. In the case of θAs and θkl , there is no correlation because generally the mass transfer rate and the solid dissolution rate are not correlated. To conclude, our Bayesian estimation approach provides quantified information about uncertainties in all the parameters as well as the relations between them. This is a big advantage of our method over the optimization which can only give a deterministic value for a parameter, even when it has a multimodal or distributed characteristic.

28


Page 28 of 62

Page 29 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5.4

Kinetics study with the proposed model parameters

We conduct sensitivity analysis using the inferred parameters (θ mean ). We perturb each de inlet sign variable at a time and see how the CO2 removal efficiency (%) = 100· n˙ inlet ˙ outlet /n˙ CO2 CO2 − n CO2 changes accordingly. As shown in Figure 11, all the design variables have significant impacts on the removal efficiency. First, when the reactor temperature increases, the CO2 removal efficiency is maintained high for longer period and the basic reactants, NaOH or Ca(OH)2, are consumed earlier. It is because a high temperature can speed up the dissolution rate of Ca(OH)2 as well as the enhancement factor E and the mass transfer rate through the rate constant k11 . Although the gas solubility, H CO2 , drops when the temperature becomes high, the increase in E would offset the decrease. However, it seems like that the removal efficiency no longer increases above a certain temperature as shown in the saturating profile in Figure 11. ˙ Second, φinlet CO2 and V , both of which are related to the inlet flow rate of CO2, have significant effects on the removal efficiency. When the inlet flow rate is low, the basic reactants, NaOH or Ca(OH)2, are consumed slowly and the CO2 removal efficiency stays high for longer period. However, the duration of high efficiency is not solely dependent on the inlet flow rate; CO2 capture rate, which is closely related to the removal efficiency, inlet ˙ ˙ decreases when φinlet CO2 or V is reduced. It is because the reduction of φCO2 and V lead to the

decrease of the driving force of the mass transfer (the concentration gradient between gas and liquid phase) as well as the overall gas holdup. However, although we can earn a high percentage of CO2 removal by maintaining low CO2 inlet flowrate, the absolute amount of captured CO2 capture would be low as well. Third, the increase of wCa(OH)2 shows an obvious effect on the removal efficiency profile. It helps to prolong the active CO2 capture period. However, in reality, it must be considered that a large loading of Ca(OH)2 could harm the mechanical operation such as impeller rotation or gas distribution. The increase of wNaOH also has similar effects unless NaOH does not exist at all. (See the dark line in the fourth plot of Figure 11). In the absence of 29



NaOH, the solution cannot completely remove the CO2 coming into the reactor. It is because the dissolution rate of Ca(OH)2 is much slower than NaOH to catch the CO2’s consumption of OH− ion. We provide one more analysis here using the obtained kinetics. In real applications, the successful operation of the process would be evaluated by how long the high removal efficiency is maintained in a single batch when fixed amount of basic materials are given and how much CO2 is captured during that time. The latter can be evaluated by how much CaCO3 is produced as well. Here two case studies, in which how they are affected by (φinlet CO2 , V˙ ) and (T , φinlet CO2 ), are done (Figure 12). In the former case, it is shown that the duration of active period increases as both design variables decrease, whereas the CaCO3 production is high at the middle of the ranges. It seems like that the southwest end region in the right upper plot of Figure 12 represents the best operating conditions. In case of (T , φinlet CO2 ), the region with dark red color in both left lower and right lower graphs, could be regarded as the best. We expect that further uncertainty analysis using uncertainty propagation techniques such as PCE, which could be one of our future studies, can help us derive more robust operating conditions.

30


Page 30 of 62

Page 31 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6

Conclusions

In this study, we developed a mathematical model for aqueous mineral carbonation process for CCU. The model describes gas-liquid mass transfer, solid dissolutions, ionic reactions and precipitation kinetics in differential algebraic equations. We introduced eight parameters in the parts where the model seems to have uncertainties. We conducted 17 sets of experiments under different conditions including one for the result validation and monitored two responses, outlet CO2 fraction and pH of the solution, in each experiment. The results from the experiments were used to estimate the posterior distributions of the parameters using Bayesian parameter estimation framework that we developed. The purpose of the Bayesian parameter estimation is to find the posterior distribution of each parameter. The posterior is defined by the product of the prior distribution of the parameter and the likelihood of data, but MCMC sampling method is necessary because the algebraic formulation of the posterior with respect to θ is unattainable. For the prior, we adopted Uniform distribution and for the likelihood, Gaussian error assumption in which iterative samplings were used to estimate the covariance vector of Gaussian error is used. As a result, we earned the stochastic model responses which fit the experimental measurements very well. It shows that our model and the parameter estimation method were capable of capturing the important features of the mineral carbonation reactor. Some of the parameters were revealed to have multimodal or distributed natures, although we designed them to have deterministic values in the model development stage. It means that a deterministic model can be extended to a stochastic model through Bayesian inference. Based on the obtained parameter estimates, we provided a rigorous analysis on the effect of each design variable on the reactor’s performance such as CO2 removal efficiency and the production of CaCO3. This can be useful information for operating the mineral carbonation process. We recommend our methodology for the parameter estimation problems which have nonideal or hard-to-predict behaviors in systems. If model errors are not fully unavoidable, it 31



is good to keep the uncertainties with the model for follwing applications. We distribute a toolbox (https://github.com/jihyunbak/BayesChemEng) with an end-to-end pipeline from creating surrogate models for dynamic responses through inferencing posterior distributions using MCMC to analyzing results with plotting tools. We believe that it can be applied in every field which requires a stochastic model to quantify uncertainties in the output response, such as robust design, stochastic optimal control, and model-based design of experiments .

32


Page 32 of 62

Page 33 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Acknowledgement This work was supported by the Energy Efficiency & Resources Core Technology Program of the Korea Institute of Energy Technology Evaluation and Planning(KETEP) granted financial resource from the Ministry of Trade, Industry & Energy(MOTIE), Republic of Korea (No. 20152010201850) and also supported by Korea Institute of Science and Technology (KIST) Institutional Program (2E29482).

Supporting Information Available The Supporting Information is available free of charge on the ACS Publications website at DOI: . • Supplementary Material • All codes are maintained in GitHub: https://github.com/jihyunbak/BayesChemEng

References (1) Bui, M. et al. Carbon capture and storage (CCS): the way forward. Energy & Environmental Science 2018, 11, 1062–1176. (2) Otto, A.; Grube, T.; Schiebahn, S.; Stolten, D. Closing the loop: captured CO2 as a feedstock in the chemical industry. Energy & Environmental Science 2015, 8, 3283– 3297. (3) Herron, J. A.; Kim, J.; Upadhye, A. A.; Huber, G. W.; Maravelias, C. T. A general framework for the assessment of solar fuel technologies. Energy & Environmental Science 2015, 8, 126–157.

33



Page 34 of 62

(4) Pan, S.-Y.; Adhikari, R.; Chen, Y.-H.; Li, P.; Chiang, P.-C. Integrated and innovative steel slag utilization for iron reclamation, green material production and CO2 fixation via accelerated carbonation. Journal of Cleaner Production 2016, 137, 617–631. (5) Sanna, A.; Uibu, M.; Caramanna, G.; Kuusik, R.; Maroto-Valer, M. M. A review of mineral carbonation technologies to sequester CO2. Chemical Society Reviews 2014, 43, 8049–8080. (6) IPCC, IPCC Special Report on Carbon Dioxide Capture and Storage, Prepared by Working Group III of the Intergovernmental Panel on Climate Change; Report, 2005. (7) Orica,

http://www.orica.com/About-Us/Innovation—Technology/mineral-

carbonation-international#.W0b2htIzY2x, accessed July 2018. (8) Futurism,

https://futurism.com/revolutionary-carbon-capture-method-makes-

building-materials-out-of-emissions/, accessed July 2018. (9) Calera Corporation, http://www.calera.com/beneficial-reuse-of-co2/scale-up.html, accessed July 2018. (10) Cleantech

Group,

https://www.cleantech.com/between-a-rock-and-hard-place-

commercializing-co2-through-mineralization/, accessed July 2018.

(11) Chosun Biz internet news, http://biz.chosun.com/site/data/html_dir/2017/06/26/2017062600894.htm accessed July 2018.

(12) The Asia Business Daily internet news, http://www.asiae.co.kr/news/view.htm?idxno=2015070916592 accessed July 2018. (13) Scott, A. Learning to love CO2. Chemical & Engineering News 2015, 93, 10–16. (14) Huntzinger, D. N.; Gierke, J. S.; Kawatra, S. K.; Eisele, T. C.; Sutter, L. L. Carbon Dioxide Sequestration in Cement Kiln Dust through Mineral Carbonation. Environmental Science & Technology 2009, 43, 1986–1992. 34


Page 35 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) Kwon, S.; Fan, M.; Dacosta, H. F. M.; Russell, A. G.; Tsouris, C. Reaction Kinetics of CO2 Carbonation with Mg-Rich Minerals. The Journal of Physical Chemistry A 2011, 115, 7638–7644. (16) Matter, J. M.; Kelemen, P. B. Permanent storage of carbon dioxide in geological reservoirs by mineral carbonation. Nature Geoscience 2009, 2, 837. (17) Mitchell, M. J.; Jensen, O. E.; Cliffe, K. A.; Maroto-Valer, M. M. A model of carbon dioxide dissolution and mineral carbonation kinetics. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 2009, (18) Montes-Hernandez, G.; Pérez-Lépez, R.; Renard, F.; Nieto, J. M.; Charlet, L. Mineral sequestration of CO2 by aqueous carbonation of coal combustion fly-ash. Journal of Hazardous Materials 2009, 161, 1347–1354. (19) Velts, O.; Uibu, M.; Kallas, J.; Kuusik, R. Waste oil shale ash as a novel source of calcium for precipitated calcium carbonate: Carbonation mechanism, modeling, and product characterization. Journal of Hazardous Materials 2011, 195, 139–146. (20) Najm, H. N. Uncertainty Quantification and Polynomial Chaos Techniques in Computational Fluid Dynamics. Annual Review of Fluid Mechanics 2009, 41, 35–52. (21) Hosder, S.; Walters, R.; Perez, R. 44th AIAA Aerospace Sciences Meeting and Exhibit; Aerospace Sciences Meetings; American Institute of Aeronautics and Astronautics, 2006. (22) Paulson, J. A.; Mesbah, A. An efficient method for stochastic optimal control with joint chance constraints for nonlinear systems. International Journal of Robust and Nonlinear Control 2017, 1–21. (23) Mesbah, A. Stochastic model predictive control with active uncertainty learning: A Survey on dual control. Annual Reviews in Control 2018, 45, 107 – 117. 35



(24) Huan, X.; Marzouk, Y. M. Simulation-based optimal Bayesian experimental design for nonlinear systems. Journal of Computational Physics 2013, 232, 288–317. (25) Kastner, C. A.; Braumann, A.; Man, P. L. W.; Mosbach, S.; Brownbridge, G. P. E.; Akroyd, J.; Kraft, M.; Himawan, C. Bayesian parameter estimation for a jet-milling model using Metropolis-Hastings and Wang-Landau sampling. Chemical Engineering Science 2013, 89, 244–257. (26) Konomi, B. A.; Karagiannis, G.; Lai, K.; Lin, G. Bayesian Treed Calibration: An Application to Carbon Capture With AX Sorbent. Journal of the American Statistical Association 2017, 112, 37–53. (27) Li, K.; Mahapatra, P.; Bhat, K. S.; Miller, D. C.; Mebane, D. S. Multi-scale modeling of an amine sorbent fluidized bed adsorber with dynamic discrepancy reduced modeling. React. Chem. Eng. 2017, 2, 550–560. (28) Coleman, M. C.; Block, D. E. Bayesian parameter estimation with informative priors for nonlinear systems. AIChE Journal 2006, 52, 651–667. (29) Bak, J. H.; Pillow, J. W. Adaptive stimulus selection for multi-alternative psychometric functions with lapses. Journal of Vision 2018, 18, 4. (30) Nelson, C.; Renate, M.; Lloyd, K.; Ben, L. Bayesian methods for cosmological parameter estimation from cosmic microwave background measurements. Classical and Quantum Gravity 2001, 18, 2677. (31) Guofang, X.; Brady, M.; Noble, J. A.; Yongyue, Z. Segmentation of ultrasound B-mode images with intensity inhomogeneity correction. IEEE Transactions on Medical Imaging 2002, 21, 48–57. (32) Ghahramani, Z. Probabilistic machine learning and artificial intelligence. Nature 2015, 521, 452. 36


Page 36 of 62

Page 37 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(33) Mosbach, S.; Braumann, A.; Man, P. L. W.; Kastner, C. A.; Brownbridge, G. P. E.; Kraft, M. Iterative improvement of Bayesian parameter estimates for an engine model by means of experimental design. Combustion and Flame 2012, 159, 1303–1313. (34) Abbasi, M.; Gholami, A. Polynomial chaos expansion for nonlinear geophysical inverse problems. Geophysics 2017, 82, R259–R268. (35) Hindmarsh, A. C.; Brown, P. N.; Grant, K. E.; Lee, S. L.; Serban, R.; Shumaker, D. E.; Woodward, C. S. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Trans. Math. Softw. 2005, 31, 363–396. (36) Kang, D.; Park, S.; Jo, H.; Park, J. Carbon fixation using calcium oxide by an aqueous approach at moderate conditions. Chemical Engineering Journal 2014, 248, 200–207. (37) Kang, D.; Lee, M.-G.; Jo, H.; Yoo, Y.; Lee, S.-Y.; Park, J. Carbon capture and utilization using industrial wastewater under ambient conditions. Chemical Engineering Journal 2017, 308, 1073–1080. (38) Brauer, H. Particle/fluid transport processes. Prog. Chem. Eng. 1981, 19, 61–99. (39) Ratcliff, G.; Holdcroft, J. Diffusivities of gases in aqueous electrolyte solutions. Trans. Inst. Chem. Eng 1963, 41, 315–319. (40) Westerterp, K. R.; Van Swaaij, W.; Beenackers, A.; Kramers, H. Chemical reactor design and operation. 1984, (41) Zhang, D.; Deen, N.; Kuipers, J. Numerical modeling of hydrodynamics, mass transfer and chemical reaction in bubble Columns. 2007, (42) Weisenberger, S.; Schumpe, A. Estimation of gas solubilities in salt solutions at temperatures from 273 K to 363 K. AIChE Journal 1996, 42, 298–300.

37



(43) Versteeg, G. F.; Van Swaaij, W. P. Solubility and diffusivity of acid gases (carbon dioxide, nitrous oxide) in aqueous alkanolamine solutions. Journal of Chemical & Engineering Data 1988, 33, 29–34. (44) Box, G. E.; Hunter, W. G.; Hunter, J. S. Statistics for experimenters; John Wiley and sons New York, 1978. (45) Box, G. E.; Draper, N. R. Empirical model-building and response surfaces.; John Wiley & Sons, 1987. (46) Atkinson, A.; Donev, A.; Tobias, R. Optimum experimental designs, with SAS ; Oxford University Press, 2007; Vol. 34. (47) Chen, J.; Wang, K.-P. Sequential Experimental Design Strategy for Optimal Batch Profiles Using Hybrid Function Approximations. Industrial & Engineering Chemistry Research 2004, 43, 5260–5274. (48) McKay, M. D.; Beckman, R. J.; Conover, W. J. Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Technometrics 1979, 21, 239–245. (49) Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A. Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics 1953, 21 . (50) Jones, D. R.; Perttunen, C. D.; Stuckman, B. E. Lipschitzian optimization without the Lipschitz constant. Journal of Optimization Theory and Applications 1993, 79, 157–181. (51) Finkel, D. E.; Kelley, C. Convergence analysis of the DIRECT algorithm. Optimization Online 2004, 14, 1–10. (52) Na, J.; Lim, Y.; Han, C. A modified DIRECT algorithm for hidden constraints in an LNG process optimization. Energy 2017, 126, 488 – 500. 38


Page 38 of 62

Page 39 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(53) Kakaraniya, S.; Gupta, A.; Mehra, A. Reactive Precipitation in Gas-Slurry Systems: The CO2-Ca(OH)2-CaCO3 System. Industrial & Engineering Chemistry Research 2007, 46, 3170–3179.

39



Figures and Tables Vent

Reactor system

Gas analyzer MFC

Chiller

CO2

Gas mixer

t

High Purity N2

pH meter

Reactor pH

Water bath

N2

Vent

Analyzer CO2 vol%

Gas preparation

Water bath

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 62

t

Reactor

Experiment design variables ( )

Experiment output ( )

Figure 1: Schematic process flow diagram of aqueous mineral carbonation process for CO2 absorption and utilization.

𝑸𝑸𝑵𝑵

𝑸𝑸𝑵𝑵−𝟏𝟏

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝑵𝑵 ]

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝑵𝑵−𝟏𝟏 ]

Well-mixed Ions and solids

𝑸𝑸𝟐𝟐 𝟏𝟏

𝑸𝑸

𝑸𝑸𝟎𝟎

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝟑𝟑 ] [𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝟐𝟐 ]

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝟏𝟏 ]

Figure 2: Schematic diagram of the compartment reactor model.

40


Page 41 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Experiment

M=16 Orthogonal Array

( 1) 2

Exp. M

Model

Design variable

Exp. 1 Parameter

2

3

time

Figure 3: Schematic diagram of the problem setting. a) Generating model

b) Iterative algorithm for posterior inference prior

Design variable

sample posterior

final posterior

Modeled response Observed response

Parameter

data (x,y) model

Error

likelihood estimate error cov.

predict response

Figure 4: a) A graphical representation of model generation. The observed response is considered as a noisy version of the modeled (ideal) response. b) Schematic for the Bayesian inference.

41



12

3

10

2

8

Page 42 of 62

1

6 0

95% CI

4

-1

2

-2

0

-3 -3

-2 0

2000 4000 6000 8000 10000 12000

-2

-1

0

1

2

3

Figure 5: a) Estimated output response (φoutlet CO2 ) fitting and b) joint marginal posterior distribution of log10 (θkl ) and log10 (θAs ) for Exp.6 using mode, mean, and optimization method.

2

2

2

2

0

0

0

0

-2

-2

-2

-2

-2

0

2

-2

0

-2

2

0

2

-2

0

2

Figure 6: Tempering control (tempering factor as 1, 4, 8, and 10 respectively) for checking multimodal posterior distribution between θAs and θkl ; also compare the θmean and θmode movements. × denotes mean and + denotes mode.

42


Page 43 of 62

10

1

log10(Error)

0.1

b) f(θf(θ

f(θopt.) f(θmean) f(θmode) fmean

opt.)

mean)

Train & Validation Set (Orthogonal array DoE)


10

Error

1

Test set

f(θmode) fmean Test set

a)

Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1


0.1 2

4

6

8

10 12 14 16

2

Experiment number

4

6

8

10 12 14 16

Experiment number

Figure 7: Error comparison through different parameter estimation methods. Visualize the error with a) CO2 volumetric concentration of outlet flue gas (φoutlet CO2 ) and b) pH

43



Figure 8: Two dimentional joint marginal posterior distribution with 20,000 final samples of MCMC. The center of dotted crosslines is θ mean .

44


Page 44 of 62

Page 45 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.1

0.1

0.1

0.05

0.05

0.05

0

0 -2

0

2

0 -2

0

2

0.1

0.1

0.1

0.05

0.05

0.05

0 0

0.1

0.2

0.3

0 -2

0.1

0.1

0.05

0.05

0 -6

-4

-2

0

0.1

0.2

0

0.5

1

0.3

0 0

2

0 -2

0

2

Figure 9: Marginal posterior distributions of each parameter from different kinds of datasets. Each gray line corresponds to the posterior which is inferred from two responses data from same experiment.

45



Page 46 of 62

-2

-2

-2

-2

-4

-4

-4

-4

-6

-6

-6

-6

-2

0

2

-2

0

-2

2

0

2

-2

0

2

Figure 10: Tempering control (tempering factor as 1, 4, 8, and 10 respectively) for checking multimodal posterior distribution between θAs and θd0s ; also compare the θmean and θmode movements. × denotes mean and + denotes mode.

100

100

100

100

100

80

80

80

80

80

60

60

60

60

60

40

40

40

40

40

20

20

20

20

20

0

0

0

0

0

2000

4000

6000

0

2000

4000

6000

0

2000

4000

6000

High

Low

0 0

2000

4000

6000

0

2000

4000

6000

Figure 11: The sensitivity analysis of design variables. While one design variable is perturbed at a time within a range of the specified low and high values, the other variables are fixed to the median value of the operating range, i.e. T to 46.5 ◦ C, φinlet CO2 to 16 %, wCa(OH)2 to 2 wt%, wNaOH to 1.5 wt%, and V˙ to 1.25 L/min. The perturbation is conducted to have same interval between the values.

46


Page 47 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Duration of active period (s) 2

104 10

1.5

2

30

1.5 5

1

20

1 10

0.5

0.5 0

10

20

0

30

10

20

30

104

30

30

30

20

20

10

10

2

20 1

10 0

40

0

60

40

60

Figure 12: The case study results using the estimated parameter values

47



Page 48 of 62

Table 1: 2-(4,5,1) orthogonal array experiment design Exp. # T (◦ C) 1 25 2 25 3 25 4 25 5 40 6 40 7 40 8 40 9 50 10 50 11 50 12 50 13 68 14 68 15 68 16 68 Test 70

φinlet wNaOH (wt%) CO2 (vol%) wCa(OH)2 (wt%) 2 1 0 10 1.5 1 20 2 2 30 3 3 2 1.5 2 10 1 3 20 3 0 30 2 1 2 2 3 10 3 2 20 1 1 30 1.5 0 2 3 1 10 2 0 20 1.5 3 30 1 2 10 3 1

V˙ (L/min) 0.5 1 1.5 2 2 1.5 1 0.5 1 0.5 2 1.5 1.5 2 0.5 1 2

Table 2: The prior range for the parameter θ Parameter log(θE ) log(θAs ) θa θb log(θc ) log(θdb ) log(θd0s ) log(θkl )

Min. -3 -3 0 0 -2 0 -7 -3

Max. 3 3 0.3 0.3 2 1.48 -1 3

Description adjustment factor for enhancement factor adjustment factor for surface area of CaCO3 related to gas holdup related to gas holdup related to gas holdup initial diameter of bubbles initial diameter of CaCO3 particles adjustment factor for mass transfer coefficient

48


Page 49 of 62

Experiment

Bayesian Inference

Kinetics Model 𝒅𝒅 𝐂𝐂𝐎𝐎𝟐𝟐 𝟐𝟐+ = 𝒇𝒇(𝐇𝐇𝐇𝐇𝐎𝐎− 𝟑𝟑 , 𝐂𝐂𝐚𝐚 , … ) 𝒅𝒅𝒅𝒅

Parameter 2

MCMC samples for posterior Posterior mode Posterior mean Optimization

CO2

Design variables Parameter 1

Hidden Parameters

Marginal distribution

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pH

Parameter 1

Parameter 2

Each experiment All experiments All experiments (tempered)

outlet 𝑃𝑃CO 2

Figure 13: Graphical abstract

49


MFC

Chiller

CO2

Gas mixer

Experiment design variables ( )

t

High Purity N2

pH meter

Reactor pH

Water bath

N2

Vent

Analyzer Gas analyzer

Water bath

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Reactor system

Page 50 of 62

CO2 vol%

Gas preparation


Vent

t

Reactor


Experiment output ( )

Page 51 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


𝑚

𝑚

Experiment

𝒙(𝑚)

M=16 Orthogonal Array

(𝒚1𝑚 , 𝒚2𝑚 ), 1 ≤ 𝑚 ≤ 𝑀

𝑚

𝒙(𝑚) = (𝑇 𝑚 , 𝜙 inlet , 𝑤Ca OH 2 , 𝑤NaOH , 𝑉ሶ (𝑚) ) CO2 𝜽 = (𝜃𝐸 , 𝜃𝐴𝑠 , 𝜃𝑎 , 𝜃𝑏 , 𝜃𝑐 , 𝜃𝑑𝑏 , 𝜃𝑑𝑠 , 𝜃𝑘𝑙 )

𝒙(1) Design variables

𝑦Ԧ

(1)

𝑦Ԧ2

(1)

Model

(1)

𝒚2

3 2

𝜽 Parameters

Exp. M

𝒚1

𝒚model

Exp. 1

𝐸(𝑓) 𝜎(𝑓)

time


Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

𝑸𝑸𝑵𝑵

𝑸𝑸𝑵𝑵−𝟏𝟏

Page 52 of 62

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝑵𝑵 ]

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝑵𝑵−𝟏𝟏 ]

Well-mixed Ions and solids

𝑸𝑸𝟐𝟐

𝑸𝑸𝟏𝟏

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝟑𝟑 ] [𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝟐𝟐 ]

[𝑪𝑪𝑶𝑶𝟐𝟐 (𝒈𝒈)𝟏𝟏 ]


𝑸𝑸𝟎𝟎

Page 53 of 62 1 2 3 4M=16 5Orthogonal 6 7Array 8 9 10 11 12 13 14 15

Industrial & Engineering Chemistry Research Experiment

( 1) 2

Design variable

Exp. 1


Parameter

Exp. M

Model

time

2

3

a) Generating model

b) Iterative algorithm for posterior inference Industrial & Engineering Chemistry Research prior

1 Design 2 variable

3 4 5 6 7 8 Parameter 9 10 11 12 13 14 15 16 17 18

Page 54 of 62

sample posterior

Modeled response Observed response

data (x,y) model

Error

likelihood estimate error cov.


predict response

final posterior

Page 55 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43


12

3

10

2

8

1

6 95% CI

4

0 -1

2

-2

0 -2 0

-3 ACS Paragon Plus Environment

2000 4000 6000 8000 10000 12000

-3

-2

-1

0

1

2

3


1

log10(Error)

0.1

b) f(θf(θ


opt.)

mean)


10

Error

1


1


0.1 2

4

6

8

10 12 14 16

Experiment number

2 ACS Paragon Plus Environment

4

6

8

10 12 14 16

Experiment number

Test set

f(θmode) fmean Test set

a)

10

Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 56 of 62

Page 57 of 62


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



0.1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.05

0

0.1

0.1

0.05

0.05

0 -2

0

2

0 -2

0

2

0.1

0.1

0.1

0.05

0.05

0.05

0 0

0.1

0.2

0.3

0 -2

0.1

0.1

0.05

0.05

0 -6

-4

-2

0

Page 58 of 62

0

0.1

0.2

0

0.5

1

0 0

ACS Paragon -2 Plus Environment 0 2

2

0.3

Page 59 of 62 1 2 3 4 5 6 7 8


-2

-2

-2

-2

-4

-4

-4

-4

-6

-6

-6

-6


-2

0

2

-2

0

2

-2

0

2

-2

0

2

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8

Page 60 of 62

2

2

2

2

0

0

0

0

-2

-2

-2

-2


-2

0

2

-2

0

2

-2

0

2

-2

0

2

Page 61 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18


100

100

100

100

100

80

80

80

80

80

60

60

60

60

60

40

40

40

40

40

20

20

20

20

20

0

0

0

0

0

0

2000

4000

6000

0

2000

4000

6000

0

2000

4000

6000


0

2000

4000

6000

High

Low

0

2000

4000

6000

Industrial & Engineering Chemistry Research Duration of active period (s) 4

10 10

2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

1.5

Page 62 of 62

2

30

1.5 5

1

20

1 10

0.5

0.5 0

10

20

0

30

10

20

30

104

30

30

30

20

20

10

10

2

20 1

10 40

0 Plus Environment ACS Paragon

60

0

40

60

Bayesian inference of aqueous mineral carbonation kinetics for

Recommend Documents