Available online at www.sciencedirect.com
ScienceDirect Energy Procedia 114 (2017) 3558 – 3563
13th International Conference on Greenhouse Gas Control Technologies, GHGT-13, 14-18 November 2016, Lausanne, Switzerland
Development of a Numerical Method to Specify the CO2 Seepage Information in the Ocean Ryosuke Sakaizawa1, Toru Sato1, Hiroyuki Oyama1, Takero Yoshida1 1
Department of Ocean Technology, Policy and Environment, The University of Tokyo 5-1-5 Kashiwanoha, Kashiwa 277-8356, Japan
Abstract Carbon Capture and Storage (CCS) is one of the important solutions against the increase of CO2 emissions in the world. Especially, in Japan, CO2 storage under the seafloor gathers attention and is promoted to reach practical stage until 2020. However, in the case of an accident, although very rare, CO2 bubbles can seep out into the seawater. For this reason, it is very convenient to have a numerical method to detect the seepage information, such as flux and location. Because measurement data detected by several sensors is used, and information of a past seepage is estimated, such a detection method is called an inverse method. Inverse methods are mainly used for atmospheric pollution research, and also for underground water or ocean research. However, these studies did not have similar conditions as the CCS accident case we predict; three dimension, continuous release and unsteady flow. Mori et al. (2014) developed the adjoint method for the CCS but this method needs further improvement. When we will conduct real estimation, there are some limitations; the number of sensors and of course cost. We need to obtain an efficient procedure such as the arrangement of measurement sensors and an effective threshold of measurement data. In this paper, we conducted the 2D test calculation and estimated the leakage information with Mori’s method. We evaluated their method and discussed the implementation. © 2017 by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license © 2017Published The Authors. Published by Elsevier Ltd. (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of GHGT-13. Peer-review under responsibility of the organizing committee of GHGT-13.
Keywords: CO2 seepage, Numerical method, Adjoint method:
1. Introduction Carbon Capture and Storage (CCS) is one of the important solutions for climate change mitigation in the world. With CCS, CO2 emitted from coal-fired power plants or other sources is captured and then is injected into an aquifer
1876-6102 © 2017 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of GHGT-13. doi:10.1016/j.egypro.2017.03.1485
Ryosuke Sakaizawa et al. / Energy Procedia 114 (2017) 3558 – 3563
3559
under the ground. This technology attracts attention because it sustains the amount of CO2 in the atmosphere. Especially in Japan, CCS under the seabed is promoted to reach practical stage until 2020 and demonstration in a coast of Japan has already started. However, there is a risk of CO2 bubbles seep into the seawater by a natural disaster such as an earthquake. In this case, we need to specify the leakage information such as the location and the amount of CO2 leakage as soon as possible to preserve the ecosystem. The detection method of leakage information is based on the data assimilation; one of the initial problems. Monte-Carlo simulation and backward methods are also data assimilation methods but they need long time and sometimes they are not stable. Moreover, in ocean modelling, because of the large scale computational domain, the above methods are not suitable. In this work, we used the adjoint method, which is one of the inverse methods. Measurement data is detected by several sensors and a past seepage information is estimated. This inverse method calculates the diffusion of adjoint probability and is able to specify the leakage information with an evaluate function. It can save more time and is more stable than other forward and backward methods. This method was also applied for specification of pollutant information in the atmosphere or inside of buildings. Neupaur (1999, 2000) [1][2] also used the adjoint method to specify the pollution source of underground water. These studies are for two dimension computational domain under steady flow condition. Liu and Zhai (2007, 2008) [3] [4] proposed an estimation method to calculate the likelihood of a gaseous pollutant source location with the adjoint method. Zhang et al. (2015) [5] succeeded to identify the location of two contaminant sources with the flux information of specific points. However, they focus on instantaneous pollution release case under steady flow. It is not like a CCS accident; CO2 seepage should be continuous and under ocean current; it is not steady flow. Zhai et al. (2012)[6] and Zhang et al. (2013)[7] studied estimation against continuous release inside structures, but there is another serious problem that needs to be addressed. This problem is leakage information; location, start time of leakage or leakage amount. Mori et al. (2014) [8] focused on the CCS leakage accident and adopted an adjoint marginal sensitivity method (AMSM) to estimate seepage information by using measurement data detected at several sensors set in a virtual oceanic field. Their method was successful to some extent. However, their method needs improvement to make an efficient procedure to detect the leakage information; the dependence on spatial arrangement of sensors, the minimum concentration for the estimation etc. Our main objective is to estimate CO2 leakage information more clearly and effectively. In this paper, we conducted some estimation cases in order to evaluate the universality; dependence on the arrangements of sensors, flow velocity etc. in a continuous CO2 release case. 2. Methodology 2.1. Theory Mori et al. (2014)[8] adopted the AMSM (Adjoint Marginal Sensitivity Method) of Dimov et al. (1996)[9] to estimate the leakage information. In short, AMSM calculates the diffusion of adjoint probability in reverse time direction with adjoint equation. In general, transport equation of CO2 can be expressed as ∂C ∂V j C ∂ ⎡ ∂ 2C ⎤ + = ⎢ν C ⎥ + M 0δ ( x − x 0 )δ (t − t 0 ) ∂t ∂x j ∂x j ⎣⎢ ∂x 2j ⎥⎦
(1)
where Vj is the velocity component in xj direction. C is the concentration of CO2 and M 0 , x0 and t 0 are the flux leakage amount, the location and time of the occurred CO2 leakage, respectively. Adjoint equation can be derived from the following forward transport equation:
3560
Ryosuke Sakaizawa et al. / Energy Procedia 114 (2017) 3558 – 3563
Fig. 1. The concept of adjoint method
∂ψ ∂V jψ ∂ − = ∂τ ∂x j ∂x j
(2)
⎡ ∂ 2ψ ⎤ + Γ0δ ( x − xi )δ (−i ) ⎢ν C 2 ⎥ ⎢⎣ ∂x j ⎥⎦
where τ is the time step in backward time, ψ is the adjoint probability and Γ0 is the adjoint probability flux (Γ0 = 1.0). xi and τ i are the location and CO2 measurement time of i th sensor. We can consider the eq. (2) as the transport equation of ψ. That means that we can estimate the leakage information by calculating the diffusion of ψ in reverse time direction (backward time direction). In the moment of CO2 release case, we consider that CO2 was released at t = t1 and measured by a sensor at t = tobs as C1 (Fig. 1). The time from CO2 leakage to measurement is defined as “travel time” and the ratio between leakage and measured amount as “marginal sensitivity”. In this example case, travel time is tobs – t1 and marginal sensitivity is C1 / M0. In backward calculation, adjoint probability was released at i th sensor and τ = τleak and calculated ψ1 at the point CO2 released in forward calculation and backward time τ = τ1. For ψ, travel time is τ1 - τleak and marginal sensitivity is ψ1 / Γ0. The main point of adjoint method is when travel time becomes equal between forward and backward, marginal sensitivity should be
C obs ψ cal C = ⇔ M 0 = obs Γ0 M 0 dt Γ0 dt ψ cal
(3)
We can consider the continuous release case as the sum of the momentum release cases and they are expressed as: For C1 (released at t1 and measured at tobs)
M0 =
For C2 (released at t2 and measured at tobs)
M0 =
C1
ψ cal (τ = τ 1 − τ leak )
Γ0
C2 Γ0 ψ cal (τ = τ 2 − τ leak )
(4)
(5)
For CN (released at tN and measured at tobs) Where
M0 =
CN Γ ψ cal (τ = τ Ν − τ leak ) 0
Cobs = C1 + C2 + + C N
From the above equations, Mori et al. (2014) derived the following equation for continuous release case;
(6)
(7)
Ryosuke Sakaizawa et al. / Energy Procedia 114 (2017) 3558 – 3563
M0 =
3561
Γ0Cobs tobs
1 / dt ∫ ψ (τ )dτ
(8)
0
2.2. Test Simulation In this study, we conducted forward calculation to simulate the diffusion of dissolved CO2 in a virtual 2D oceanic site and used the results as measurement data for backward calculation. This is because we can’t conduct CO2 seeping examination in the ocean due to regulations. In the forward calculation, our calculation was based on the Marine Environmental Committee (MEC) Ocean Model (e.g. Sato et al., 2006; Kano et al., 2010) [9] [10]. The governing equations are the Navier-Stokes equations using hydrostatic approximation and continuous equation. However, in this study, we used the following flow field to express tidal flow in the test computational domain.
u = 0.5 × cos(2πt / 44714.2)
(9) (10)
v = 0.5 × sin(π + 2πt / 44714.2)
where u, v are the velocity components of the tidal flow in x, y directions. The CO2 transportation is governed by the following equation; (11)
∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ DC ⎟ + ⎜ AV = − ⎜ AH ⎟ + ⎜ AH ⎟ Dt ∂x ⎝ ∂x ⎠ ∂y ⎜⎝ ∂y ⎟⎠ ∂z ⎝ ∂z ⎠
The 2D computational domain of the test calculation has 60×60 grids in horizontal direction and each grid size is 1000 m × 100 0m in the horizontal plane and 20 m in the vertical direction (Fig. 2). The forward calculation conducted for 200,000 time steps equivalent to 200,000 sec. CO2 seepage started at 100,000 seconds later at (30, 30) and the amount was 1.0 kg/m3/sec continuously. We arranged 24 sensors in the computational domain and measurement and estimation are conducted by each sensor. They calculated M0(x, τ) assuming each point and time step as leakage source and start time. We evaluate the error between the estimations with an evaluate function F(x, τ) and the x and τ caused the smallest F (i.e. the most acceptable estimation).
Red dot: Leakage point, Blue dot: CO2 sensors Fig. 2. The concept of adjoint method
Fig. 3. Time-series of CO2 concentration at sensor 3
3562
Ryosuke Sakaizawa et al. / Energy Procedia 114 (2017) 3558 – 3563 Table 1. Measurement data at each sensor in continuous release case tobs Cobs tobs No. (x, y) No. (x, y) 3 [step]
1 2 3 4 5 6 7 8
(20, 40) (25, 40) (30, 40) (35, 40) (40, 40) (40, 35) (40, 30) (40, 25)
172905 175787 178884 182653 185440 187000 190042 193796
[kg/m ]
1.77E+00 8.38 E+01 2.95 E+02 8.39 E+01 3.42 E -01 7.36 E+01 2.87 E+02 8.23 E+01
Cobs
No.
5.04 E -01 1.43 E+02 2.50 E+01 1.09 E-01 2.69 E+01 2.03 E+02 6.47 E+01 6.47 E+01
17 18 19 20 21 22 23 24
[kg/m3]
[sec]
9 10 11 12 13 14 15 16
(40, 20) (35, 20) (30, 20) (25, 20) (20, 20) (20, 25) (20, 30) (20, 35)
196689 198222 156514 160539 157530 164394 167364 171051
Fig.4 The time-series data of evaluate function F(x, τ)
F ( x, τ ) =
( M 0i ( x,τ ) − M 0 j ( x,τ )) 2 /(M 0i ( x,τ ) + M 0 j ( x,τ )) 2
(x, y) (25, 35) (30, 35) (35, 35) (35, 30) (35, 25) (30, 25) (25, 25) (25, 30)
tobs
[sec]
172068 174195 184195 185533 195450 196908 160852 162738
Cobs
[kg/m3]
9.27 E+02 1.19 E+03 9.96 E+02 1.20 E+03 1.00 E+03 1.26 E+03 6.52 E+02 1.09 E+03
Fig.5 Time-series of CO 2 and ψ dτ
(12)
n(n − 1) / 2
3. Results Table 1 shows the results of the forward calculation (simulated CO2 diffusion) for the coordinates of each sensor and measurement data. In Fig. 3 the time-series of CO2 concentration at sensor 3 are presented. We used the highest concentration during forward calculation as Cobs. The sensors at arranged corner; sensor 1, 5, 9, 13 measured smaller concentration compared to the other sensors and the sensors closed to leakage source (sensor 17 – 24) measured higher concentration. We used these results for backward calculation and estimated the leakage information. The time-series data of the evaluation function F(x, τ) is shown in Fig. 4. With surrounding CO2 leakage start time τ = 100,000 we got the smallest F(x, τ) at (30, 30). We can say that the leakage location is estimated constantly. The smallest F(x, τ) is obtained at τ = 101153 and (30, 30). However, in this condition the amount of CO2 leakage flux is estimated to be 6.098E-01 kg/m3s as the average of 24 sensors. When we gave 1.0E+00 kg/m3s leakage amount in the forward calculation and this estimation had 39.0% error. Fig. 5 shows the time-series of CO2 concentration and ψ dτ at sensor 3 with adjusting time direction and leakage start time. With the concept of adjoint method, CO2 travel time is 78884sec (the green line on Fig. 5) and ψ dτ should be equal to CO2 concentration. However, there is a difference in the amount because we modified the negative value of CO2 and ψ in the calculation. On Fig. 5, ψ dτ closed to CO2 concentration at 34198 sec and sensor 3 estimated the amount of CO2 flux to be 9.64E-01 (3.6% error). However, at sensor 1, ψ dτ became higher than the CO2 concentration and the estimation error of leakage flux is at least 73.4%.
Ryosuke Sakaizawa et al. / Energy Procedia 114 (2017) 3558 – 3563
3563
4. Conclusion In this paper, we used Mori’s adjoint method in a 2D test calculation model and estimated the leakage rate and location. We can estimate the location and start time of leakage accurately but the amount of leakage flux is not enough. If the factor gets a negative value in the transfer calculation, the error becomes bigger. However, depending on the measurement time tobs, this estimation can give accurate values. In future work, we need to comprehend how to choose measurement data used for estimation or implement the adjoint method directly.
References [1] R. M. Neupauer, J. L. Wilson, “Adjoint method for obtaining backward-in-time location and travel time probabilities of conservative ground water contaminant.” Water Resources Research, vol. 35, pp. 3389-3398, 1999. [2] R. M. Neupauer, B. Borchers, J. L. Wilson, “Comparison of inverse methods for reconstructing the release history of a groundwater contamination source.” Water Resources Reserch, vol. 36, pp. 2469-2475, 2000. [3] X. Liu, and Z. Zhai, “Inverse modeling methods for indoor air pollutant tracking literature review and fundamentals” Indoor Air, vol. 17, pp. 419-438, 2007. [4] X. Liu, and Z. Zhai, “Location identification for indoor instantaneous point contaminant source by probability-based inverse Computational Fluid Dynamics modeling,” Indoor Air, vol. 18, pp. 2-11, 2008. [5] T. Zhang, H. Zhou, and S. Wang, “Inverse identification of the release location, temporal rates, and sensor alarming time of airborne pollutant source,” Indoor Air, vol 25, pp. 415-427, 2015 [6] Z. Zhai, X. Liu, and H. Wang, Y. Li, J. Liu, “Experimental verification of tracking algorithm for dynamically-releasing single indoor contaminant,” BUILD SIMUL, vol 5, pp 5-14, 2012 [7] T. Zhang, S. Yin, and S. Wang, “An inverse method based on CFD to quantify the temporal release rate of a continuously released pollutant source,” Atmospheric Environment, vol. 77, pp. 62-77, 2013. [8] C. Mori, T. Sato, Y. Kano, H. Oyama, J. Kita, D. Tsumune, and Y. Maeda, “Numerical Method for Detecting the CO2 Leakage Source in the Ocean,” Energy Procedia, vol. 63, pp. 3746-3750, 2014 [9] I. Dimov, U. Jaekel, and H. Vereecken, “A Numerical Approach for Determination of Sources in Transport Equations,” Computers and Mathematics with Applications, vol. 32, pp. 31-42, 1996. [10] T. Sato, K. Tonoki, T. Yoshikawa, and Y.Tsuchiya, “Numerical and hydraulic simulations of the effect of Density Current Generator in a semi-enclosed tidal bay,” Coastal Engineering, vol. 53, pp. 49-64, 2006. [11] Y. Kano, T. Sato, J. Kita, S. Hirabayasi, S. Tabeta, “Multi-scale modeling of CO2 dispersion leaked from seafloor off the Japanese coast,” Marine Pollution Bulletin, vol. 60, pp. 215-224, 2010.