- Email: [email protected]
Sensitivity analysis of recovery efficiency in high-temperature aquifer thermal energy storage with single well
Energy 90 (2015) 1349e1359
Contents lists available at ScienceDirect
Energy journal homepage: www.elsevier.com/locate/energy
Sensitivity analysis of recovery efficiency in high-temperature aquifer thermal energy storage with single well Jun-Seo Jeon a, 1, Seung-Rae Lee a, *, Lisa Pasquinelli b, Ida Lykke Fabricius b a b
Department of Civil and Environmental Engineering, KAIST, Daejeon 305-701, South Korea Department of Civil Engineering, DTU BYG, Denmark
a r t i c l e i n f o
a b s t r a c t
Article history: Received 21 November 2014 Received in revised form 4 May 2015 Accepted 19 June 2015 Available online 17 July 2015
High-temperature aquifer thermal energy storage system usually shows higher performance than other borehole thermal energy storage systems. Although there is a limitation in the widespread use of the HTATES system because of several technical problems such as clogging, corrosion, etc., it is getting more attention as these issues are gradually alleviated. In this study, a sensitivity analysis of recovery efficiency in two cases of HT-ATES system with a single well is conducted to select key parameters. For a fractional factorial design used to choose input parameters with uniformity, the optimal Latin hypercube sampling with an enhanced stochastic evolutionary algorithm is considered. Then, the recovery efficiency is obtained using a computer model developed by COMSOL Multiphysics. With input and output variables, the surrogate modeling technique, namely the Gaussian-Kriging method with Smoothly Clopped Absolute Deviation Penalty, is utilized. Finally, the sensitivity analysis is performed based on the variation decomposition. According to the result of sensitivity analysis, the most important input variables are selected and confirmed to consider the interaction effects for each case and it is confirmed that key parameters vary with the experiment domain of hydraulic and thermal properties as well as the number of input variables. © 2015 Elsevier Ltd. All rights reserved.
Keywords: High-temperature aquifer thermal energy storage Sensitivity analysis Gaussian Kriging method Computational experiment Gassum formation
1. Introduction Energy storage technology is an essential technology for efficient energy management and utilization. Typical energy storage is to transform electric energy obtained by fossil fuel, nuclear power, and renewable source to mechanical, chemical energy form and then save it. It has advantages to reduce energy use, greenhouse gas emission, and cost for heating and cooling. Thermal energy storage is the storage of energy in the form of heat, and hence it may be directly utilized for industrial heat, as well as providing heating or power generation. The thermal energy storage functions to compensate for the temporal and quantitative differences in the supply and demand of energy, and energy saving through it and to improve the energy use efficiency. UTES (Underground thermal energy storage) systems are generally divided into ATES (aquifer thermal energy storage), CTES (cavern thermal energy storage), and BTES (borehole thermal
* Corresponding author. Tel.: þ82 42 350 3617; fax: þ82 42 350 7200. E-mail address: [email protected] (S.-R. Lee). 1 Current address: Department of Civil and Environmental Engineering, KAIST, 291, Gwahakro, Yuseong-gu, Daejeon 305-701, South Korea. http://dx.doi.org/10.1016/j.energy.2015.06.079 0360-5442/© 2015 Elsevier Ltd. All rights reserved.
energy storage). In ATES systems, thermal energy is stored in the ground water of an aquifer and groundwater is injected and extracted by one or more wells. CTES systems utilize underground cavern to store thermal energy using hot water and therefore thermal stratification has a key role in performance of the system. BTES systems are the most commonly used and thermal energy is stored in the storage medium, usually bedrock, using closed heat exchange pipe. The ATES system has higher system performance than the BTES system and any other systems using low temperature geothermal heat because of directly using groundwater with relatively high volumetric heat capacity [1]. Andersson et al. [2] also showed similar results through investigation on the performance factor, energy saving and payback time for ATES systems in Sweden (Table 1). In this regard, ATES systems are more suitable for large scale systems than BTES systems and mainly are used for seasonal thermal energy storage both of heating and cooling. Many applications of ATES are employed in European countries such as Netherland, Sweden and Germany. In particular, five systems in 1990's, 214 systems in 2000's, and 1300 ATES systems in 2010's so far have been registered in Netherlands. The Netherlands has a leading position in the world with rapid growth of installing ATES system [3] (Fig. 1).
1350
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
Table 1 Benefits in terms of economic and energy saving (Andersson et al. [2])*. System application
Performance factor
Energy saving (%)
Payback (year)
Direct heating and cooling HP supported heating and cooling HP supported heating only Direct cooling only
20e40
90e95
0e2
5e7
80e87
1e3
3e4
60e75
4e8
20e60
90e97
0e2
*Direct heat and cooling could be considered as HT-ATES and second and third system could be considered as LT-ATES.
ATES is generally divided into HT-ATES (high temperature ATES) system and LT-ATES (low temperature ATES) system. In case of HTATES system, the high temperature can be used directly for heating. It helps to minimize or eliminate the usage of heat pump, and then more energy is saving than LT-ATES system, which generally use heat pumps to increase the temperature. Despite the significant advantages, the HT-ATES is hard to employ in practice due to severe operating problems. It covers the clogging due to particle, gas bubbles and precipitation of minerals, and corrosion of components in the ground water system, etc [4]. According to the definition of high temperature (minimum storage loading temperatures on the order of 50 C) in ECES Annex 12, there are only two operational HT-ATES systems in Germany. These are installed in Neubrandenburg by 2004 with 80 C storage loading temperature [5] and at the Reichstag building in Berlin by 1999 with 70 C storage loading temperature [6]. However as mentioned issues are gradually alleviated, the application of HT-ATES is increasing especially in Germany. To predict and understand the behavior of sophisticated problem, numerical simulation is usually used. The simulator also can be used to assess the feasibility of the project, as well as to help in storage design and optimization of well locations, monitoring, and well steering [7]. Based on numerical simulation, sensitivity analysis is one of the main interests in decision-making processes (e.g. storage design, optimization of well locations) because it makes possible to quantify the leverage of each input variable to the output [8]. Therefore, sensitivity analysis is widely employed in the many engineering fields where dealing the complex problem
related to energy, geology, and geomechanics. For instance, sensitivity analysis is used for gas production from Class I hydrate reservoir by depressurization [9], robust design building energy systems [10], assessment of energy-related CO2 emissions [11], subsoil parameter estimation in mechanized tunneling [12], volcanic source modeling quality assessment and model selection [13], and sedimentary basin geothermal system [14]. Several studies on sensitivity analysis of ATES were also recently introduced. Schout et al. [15] suggested a method to estimate the recovery efficiency based on the Rayleigh number using results of sensitivity analysis. According to their results, performance of HT-ATES, which can be calculated in terms of recovery efficiency, depends on the hydraulic and operation properties (eg. permeability, injection temperature and injection volume). Bridger and Allen [16] investigated the influence of geologic layering on heat transport and storage in ATES. Results of the sensitivity analysis state that vertical anisotropy and hydraulic gradient are most sensitive on heat transport in a nonlayered aquifer, which is analogous to our case. Kim et al. [1] investigated the effect of injection rate, hydraulic conductivity, and distance between injection and production wells, and then they showed that system performance depends on the distance between boreholes, hydraulic conductivity of an aquifer, and the production/injection rate. Yapparova et al. [7] also performed numerical simulations to assess the effect of groundwater flow and injection temperature. Despite the results from previous several studies, the most sensitivity factor is different from ours. Even though the sensitivity analysis was used to select important factors, they did not confirm the interaction effects. Furthermore the onceat-time method was used in fractional factorial design corresponding to maximum and minimum values of each factor. It has the possibility to diminish the uniformity for sampling. In our study, we investigated main and interaction effects on the efficiency recovery of HT-ATES with single well using a proper computer experiment procedure. We considered two different cases and variations in hydraulic, thermal, and operation properties. Our findings will be helpful for the preliminary design and operation phases of the HT-ATES (high temperature-aquifer thermal energy system). This paper is structured with a flowchart illustrated in Fig. 2. First, the development of computational model for HT-ATES with single well is conducted using commercial program. Through the comparison between the recovery efficiency and tilt angle in thermal contours obtained by field experimental data and numerical simulation data, the developed model is validated.
Fig. 1. ATES systems in the Netherlands (Godschalk and Bakerma [3]).
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1351
Fig. 2. Flow chart for sensitivity analysis of HT-ATES.
Second, the fractional factorial sampling is completed based on Latin Hypercube sampling with optimization and criterion. Using that, the recovery efficiency for HT-ATES is obtained by the developed model. Third, the Gaussian Kriging model with penalized likelihood and SCAD penalty is employed to construct the metamodel. In the last, the sensitivity analysis using the Sobol's indices is performed and conclusions are presented.
2. Model construction and validation The commercial program COMSOL Multiphysics was used for modeling of the HT-ATES. The program provides a variety of physics interfaces which can be combined to handle multiphysics applications. For modeling of the HT-ATES, the fluid flow and heat transfer equations are considered as governing equations. The twodimensional axial symmetric model was considered by employing heat transfer module and Darcy's law in porous media implemented in COMSOL. The hydraulic gradient of HT-ATES is assumed to be zero since the aquifer is located below 1.5 km from the ground surface. In order to limit the number of input variables and see the change in physical properties, single injection/production well was considered. Physical and thermal properties of fluid and aquifer for the modeling are listed in Table 2. Input properties were chosen equal to those of the base case that was used in the work of Buscheck [17] for model validation later. Fluid density, viscosity, and specific heat variations resulting from the high injection water temperature into the aquifer were accounted for the simulations (Fig. 3). After studying the influence of different mesh sizes and considering the computation time, a discretized mesh with a total number of 9438 triangular elements, was adopted (Fig. 4). The main aspect controlling the feasibility of HT-ATES is the recovery efficiency. The recovery efficiency is defined as
ε¼
Vp $Cw $ T p Ta T p Ta ¼ Ti Ta Vi $Cw $ðTi Ta Þ
(1)
where V represents the amounts of water, Cw is the volumetric heat capacity of water, T is the average temperature, T is the temperature, and subscripts p, i, a represent the production, injection, and ambient. According to Gutierrez-Neri et al. [18], the recovery efficiencies were calculated based on the produced water temperature at
Table 2 Input parameter values used in the base case. Parameter
Units
Domain Aquifer
Aquitard
Horizontal permeability Vertical permeability Porosity Thermal conductivity Specific heat Ambient temperature Thickness
mD mD e W/(m K) J/(kg K) C m
15,000 1500 0.25 2.29 696 20 21
0.15 0.015 0.35 2.56 696 20 9
Parameter
Units
Operation Injection/Production
Storage and rest
Pumped volume Injection flow rate Injection temperature Duration
m3 m/s C days
55,000 5.371 105 90 90/90
e e e 90/90
fourth year operation for stabilization. Assuming 360 days in a year, 90 days of injection, 90 days of storage, 90 days of production, and 90 days of storage was considered as one cycle. Then the analysis was performed for four cycles. To validate the model, the results of Buscheck [17] and Schout et al. [15] were used. Schout et al. [15] also verified their model by comparing the angle of tilt for 14 scenarios and recovery efficiencies with those of Buscheck [17]. Mean thermal contour after first injection complete (90 days) and storage (180 days) for AL55 and AI90, which is computed by the developed model is shown in Fig. 5. The dotted line represents the angle of tilt, which is the angle of the mean temperature contour with the vertical. The average difference in the angles of tilt obtained by the developed model from those of Buscheck [17] and Schout et al. [15] are 3.45% and 2.80% (Table 3). The developed model results are in good agreement with the results found by Buscheck [17] and Schout et al. [15]. Therefore, the developed model for the ATES system using COMSOL could give a reasonable prediction of the recovery efficiencies.
3. Designing and modeling for computer experiments A general computer experiments process usually covers choosing factors and experimental domain, constructing a design of
1352
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
Fig. 3. Input parameter variation with temperature.
Fig. 4. Two-dimensional axial symmetric ATES model.
experiment matrix and metamodel searching, and interpreting the model and performing a sensitivity analysis (Fang et al. [8]). Based on the computer model, experimental domain is designated to seek which input factors contribute most to a response variable, here the recovery efficiency. In this study, two types of case were considered with different locations and input parameters. In Case 1, all input parameters and geometry are the same as those of Buscheck [17]. In Case 2, the sandstone layer of Gassum formation is considered as a potential layer for the application of ATES system in Denmark. The core specimens were used in the laboratory testing for collecting the input parameters. 8 and 10 input parameters were selected for Case 1 and Case 2, respectively, which may affect the output and the range was set from the minimum to the maximum value (Table 4). 3.1. Design of experiments domain All the level-combinations of the factors provide the most reliable design to estimate the main effects and some interaction
effects among factors. However, the number of runs for all the levelcombinations increases exponentially with the number of factors. Therefore, the fractional factorial design which is a subset of all the level-combinations is usually adopted. The selected subset should have a good representation of the whole combinations. Optimization rules listed in Table 5 have been widely used in searching for space-filling designs. Among them, the local search and threshold accepting have possibility to reach and stay in the local minimum. To address this problem, Kirkpatrick et al. [19] proposed the SA (simulated annealing) algorithm for a cost function that may possess many local minima. Jin et al. [20] also suggested the ESE (enhanced stochastic evolutionary) algorithm. Since ESE algorithm uses the two parameters which control the acceptance, it can adjust itself to suit different experimental design algorithms. SA and ESE algorithms are thus adopted to reach the global minimum in the optimization. Optimality criteria are also adopted to determine which subset is the best between the old and new designs. Fang et al. [21] proved that both star discrepancy and
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1353
Fig. 5. Mean thermal contour after first injection complete (90 days) and storage (180 days) for AL55 and AI90. With (a) kha ¼ 15D, kva ¼ 1:5D, khc ¼ 1:5D, kvc ¼ 0:15D, lc ¼ 2.56 W/ (m K), and injection temperature of 55 C, (b) kha ¼ 15D, kva ¼ 1:5D, khc ¼ 15,105 D, kvc ¼ 1:5,105 D, lc ¼ 0 W/(m K), and injection temperature of 95 C, superscripts h, v represents horizontal and vertical, subscript a, c indicates aquifer and aquitard.
Table 3 Angles of tilt obtained by Buscheck [17], Schout et al. [15] and COMSOL model. Scenario 90 Days
AI90 AT90 AL90 AI55 AT55 AL55 BI90 BT90 BL90 BI55 BT55 BL55 CI55 CT55 CL55
star L2 discrepancy are unsuitable for searching the best subset. Hickernell [22,23] suggested the centered discrepancy and wraparound discrepancy which satisfy the requirements of uniformity. Centered and wrap-around discrepancy are calculated as
180 Days
aB
aS
amodel DaBm DaSm aB
aS
amodel DaBm DaSm
46.1 40.3 36 16.6 13.9 13 71 68 66.4 44 39.2 35 71.3 67.1 67
49 43 38 13 14 14 75.5 75 73.5 47 38 33 68 66 62
48.01 1.91 45 4.7 40.31 4.31 19.02 2.42 17 3.1 14.03 1.03 75.42 4.42 73 5 68.19 1.79 47.23 3.23 43.6 4.4 48 13 73.3 2 71.35 4.25 68.19 1.19
65.5 65.5 64.5 34 34 29 84.5 85 85 66.5 66 64.5 84 84 84
66.37 70 66.25 36.52 31.7 27.75 84.28 82.6 84.2 67.7 68.19 65.77 85.71 84.57 82.87
0.99 2 2.31 6.02 3 0.03 0.08 2 5.31 0.23 5.6 15 5.3 5.35 6.19
64.6 65.5 63.9 32.2 23 24 82.8 80.9 81.7 64.5 64.9 64.1 83.7 80.8 81
1.77 4.5 2.35 4.32 8.7 3.75 1.48 1.7 2.5 3.2 3.29 1.67 2.01 3.77 1.87
0.87 4.5 1.75 2.52 2.3 1.25 0.22 2.4 0.8 1.2 2.19 1.27 1.71 0.57 1.13
ðCDðDn ÞÞ2 ¼
2 1 2 n Y s 13 2X 1 1 þ xkj 0:5 xkj 0:5 12 n 2 2 k¼1 j¼1 n X n Y s 1 X 1 1 þ 2 1 þ jxki 0:5j þ xji 0:5 2 2 n k¼1 j¼1 i¼1 1 xki xji 2 (2)
ðWDðDn ÞÞ2 ¼
s n Y s 4 1 X 3 xki xji 1 xki xji þ 2 3 2 n i¼1 k;j¼1
(3)
1354
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359 Table 6 Best subset of input variables.
Table 4 Factors and experimental domain for Case 1 and Case 2.
Case 1
Case 1 Parameters
Units
Domain
N
X1
X2
X3
X4
X5
X6
X7
X8
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.325 0.575 0.425 0.775 0.975 0.625 0.525 0.475 0.025 0.125 0.925 0.675 0.275 0.175 0.725 0.825 0.225 0.075 0.875 0.375
0.325 0.425 0.025 0.625 0.575 0.175 0.975 0.725 0.675 0.075 0.825 0.775 0.875 0.525 0.925 0.275 0.475 0.225 0.125 0.375
0.825 0.075 0.275 0.425 0.925 0.675 0.475 0.975 0.625 0.875 0.125 0.325 0.375 0.025 0.775 0.175 0.725 0.225 0.575 0.525
0.975 0.875 0.125 0.325 0.825 0.625 0.775 0.175 0.275 0.425 0.375 0.925 0.475 0.025 0.075 0.575 0.675 0.725 0.525 0.225
0.575 0.975 0.625 0.025 0.675 0.525 0.225 0.425 0.775 0.175 0.475 0.325 0.925 0.275 0.825 0.125 0.075 0.725 0.875 0.375
0.825 0.075 0.425 0.375 0.475 0.225 0.575 0.975 0.625 0.125 0.175 0.675 0.875 0.725 0.275 0.925 0.525 0.325 0.775 0.025
0.425 0.575 0.775 0.375 0.725 0.925 0.975 0.275 0.875 0.625 0.325 0.075 0.675 0.475 0.525 0.825 0.175 0.225 0.125 0.025
0.825 0.625 0.925 0.025 0.225 0.075 0.775 0.575 0.675 0.475 0.875 0.375 0.125 0.275 0.425 0.525 0.975 0.325 0.725 0.175
Aquifer
Aquitard
Horizontal permeability Anisotropy Porosity Thermal conductivity Specific heat Thickness Ambient temperature
mD e e W/(m K) J/(kg K) m C
[3000e75,000] [1e10] 0.25 [1.145e4.58] [600e1200] [10.5e42] 20
[7.5e–4~4.5e-3] 10 0.35 2.56 696 9 20
Parameters
Units
Operation Injection/Production
Storage and rest
m m/s C
[18,300e165,000] [1.65e-5e1.87e-4] [50e90]
e e e
Parameters
Units
Domain Aquifer
Aquitard
Horizontal permeability Anisotropy Porosity Thermal conductivity Specific heat Thickness Ambient temperature
mD e e W/(m K) J/(kg K) m C
[100e3000] [1e10] 0.25 [1.0e5.0] [600e1200] 21 55
[7.5e-4e4.5e-3] [1e10] 0.25 [1.0e5.0] [600e1200] 30 55
Parameters
Units
Operation
Pumped volume Injection flow rate Injection temperature
m3 m/s C
Pumped volume Injection flow rate Injection temperature
3
Case 2
Injection/Production
Storage and rest
[10,000e100,000] [1e-5e10e-5] [65e95]
e e e
Table 5 Considered optimization and criterion for Latin Hypercube sampling. Optimization
Criterion Distance
Uniformity
Threshold Accepting
Minimax/ maximin fp -criterion e e
Star discrepancy
Local Search Simulated Annealing Enhanced Stochastic Evolutionary
Star L2 discrepancy Centered discrepancy Wrap-around discrepancy
The discrepancy for the best subset selected by various optimality criteria are listed in Table 5. With any optimality, the ESE algorithm provides better uniformity than the SA algorithm. Thus the best subset is determined by the ESE with the centered discrepancy and wrap-around discrepancy. The design of experiments with n ¼ 20 and s ¼ 8 for Case 1, and n ¼ 20 and s ¼ 10 for Case 2 are listed in Table 6. The selected factors are uniformly scattered as shown in the scatter plot of LHS (20, 8) and LHS (20, 10) (Fig. 6). 3.2. Modeling of metamodel A computer model is an alternative of field and laboratory testing to simulate the complicated physical phenomena and to study any relationships among various parameters. However, the computational efforts increase dramatically for design optimization and probabilistic analysis based on the computer model. To alleviate these problems, the simpler and much faster to compute than the true model with accuracy is adopted by many researchers [24e26]. Kleijnen [27] called such an approximate model “model of
Case 2 N
X1
X2
X3
X4
X5
X6
X7
X8
X9
X10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.275 0.075 0.675 0.225 0.525 0.625 0.925 0.175 0.325 0.375 0.425 0.975 0.025 0.725 0.875 0.825 0.575 0.475 0.125 0.775
0.825 0.175 0.075 0.225 0.025 0.125 0.325 0.525 0.375 0.875 0.575 0.975 0.625 0.725 0.425 0.675 0.475 0.775 0.925 0.275
0.375 0.175 0.775 0.975 0.275 0.525 0.625 0.075 0.675 0.425 0.925 0.125 0.575 0.225 0.325 0.825 0.025 0.875 0.725 0.475
0.075 0.625 0.575 0.225 0.325 0.775 0.025 0.875 0.125 0.825 0.425 0.675 0.375 0.175 0.275 0.725 0.525 0.975 0.475 0.925
0.925 0.525 0.675 0.625 0.075 0.975 0.475 0.125 0.175 0.425 0.375 0.575 0.775 0.225 0.725 0.025 0.875 0.825 0.325 0.275
0.625 0.375 0.825 0.425 0.325 0.475 0.975 0.775 0.175 0.875 0.675 0.575 0.925 0.725 0.025 0.525 0.125 0.275 0.075 0.225
0.425 0.325 0.075 0.525 0.775 0.975 0.875 0.825 0.125 0.175 0.675 0.475 0.275 0.625 0.575 0.225 0.025 0.725 0.925 0.375
0.325 0.375 0.625 0.025 0.125 0.825 0.275 0.675 0.925 0.075 0.875 0.975 0.775 0.575 0.725 0.175 0.225 0.525 0.425 0.475
0.875 0.025 0.775 0.375 0.825 0.275 0.575 0.425 0.175 0.225 0.975 0.725 0.475 0.075 0.325 0.525 0.625 0.125 0.675 0.925
0.775 0.925 0.675 0.175 0.475 0.275 0.575 0.725 0.525 0.375 0.975 0.425 0.075 0.125 0.875 0.825 0.225 0.625 0.325 0.025
the model” or “metamodel”. Since the result of computer experiment is the determinant which does not contain any random errors, the metamodel is viewed as the regression on data with no random errors. In statistics, the regression method covers from the simple linear least squares method to the sophisticated neural network method. Simpson et al. [28] presented a survey by comparing various metamodel methods such as polynomial models, neural network method, and Kriging model. They recommended that the Kriging method is the best for a deterministic model, yet highly nonlinear with a moderate number of factors. Since most of computer models are defined based on governing equations dealing with the complex physical phenomena, they show usually nonlinear responses. Therefore, the Kriging method was selected for the construction of a metamodel in this study. After introducing the Kriging method by Krige [29], many researchers further developed his work. Matheron [30] suggested the Gaussian Kriging method for modeling spatial data in geostatistics. The Gaussian Kriging model is defined as
yðxÞ ¼
L X j¼0
bj Bj ðxÞ þ zðxÞ
(4)
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1355
Fig. 6. Scatter plot of the best subset.
where {Bj(x), j ¼ 1,…,L} is a basis over the experimental domain, and z(x) is assumed to be a Gaussian process with zero mean, variance and correlation function as given by
( rðq; s; tÞ ¼ corrðzðsÞ; zðtÞÞ ¼ exp
n X
) qk jsk tk j2
(5)
k¼1
where qk is the unknown smoothing parameter, n is the number of design variables, and sk, tk are kth kth components of the samples. With the assumption in the Kriging theory, the response function is composed of a regression model and stochastic process. Therefore, the first term of R.H.S in Eq. (4) represents the regression model which yields a global approximation of the design space and the second term of R.H.S in Eq. (4) represents the stochastic process which creates localized deviations. The predicted estimates at the untired points x can be written as
b b y ðxÞ ¼ bðxÞ b b þ r 0 ðxÞR1 ðqÞ y B b
(6) 0
where R is the correlation matrix, and r(x)¼(r(q;x1,x),,,,,r(q;xn,x)) . By maximizing the MLE (Maximum Likelihood Estimator), optimum unknown parameters, q and b could be obtained. The procedure of maximizing MLE are as follows (Fang et al. [8]): Step 1. Set the initial value of b to be (B0 B)1By, the least square estimator of b Step 2. Update (s2,q) by using
0
b b R1 ðqÞ y B b s 2 ¼ n1 y B b b
s
X Q b; s2 ; q ¼ [ b; s2 ; q þ n pl qj
where pl($) is a penalty function with a regularization parameter l. SCAD (Smoothly clopped absolute deviation) penalty was proposed by Fan and Li [32].
ðal qÞþ Iðq > lÞ for some a > 2 and q > 0 p0l ðqÞ ¼ l Iðq lÞ þ ða 1Þl (12) Even though l and a are unknown, they suggested that a ¼ 3.7 resulted in a good practical performance. Therefore, the Gaussian Kriging model based on the penalized likelihood with SCAD penalty was used in this study (Table 7). 4. Sensitivity analysis based on variation decomposition The metamodel g(x) can be represented in a functional ANOVA form.
gðxÞ ¼ g0 þ
(10)
X X gi ðxi Þ þ gij xi ; xj þ ,,, þ g1,,,s ðx1 ; …; xs Þ i
R1 0
i
gi1 ;…;it ðxi1 ; …; xit Þdxk ¼ 0, k ¼ i1,…,it.
(13)
Table 7 Parameters of Gaussian Kriging model with SCAD penalty. Parameters
(8)
(9)
(11)
j¼1
where
Step 3. Update b using
1 b B0 R1 ðqÞy b mle ¼ B0 R1 ðqÞB
Li and Sudjianto [31] proposed the penalized likelihood approach to alleviate a serious problem which occurred near the optimum value of parameter based on the MLE. The penalized likelihood is defined as
(7)
Then solve this equation by a numerical iteration algorithm such as NewtoneRaphson algorithm.
. v[ b; s2 ; q vq ¼ 0 1
n [ b; s2 ; q ¼ log s2 logjRðqÞj 2 2
0
1 b 2 y B b R1 ðqÞ y B b b 2s
Step 4. Iterate Step 2 and Step 3 until it converges
b l m b b s2 c q1 c q2 c q3 c q4 c q5 c q6 c q7 c q8 c q9 d q10
SCAD penalty Case 1
Case 2
0.44 0.614 0.026 8.42E-07 4.81E-06 5.65E-07 3.87E-04 2.93E-06 1.32E-05 3.07E-05 2.67E-23 e e
0.295 0.7663 0.0181 9.92E-06 3.20E-03 7.671 3.20E-03 7.20E-03 7.20E-03 4.79E-05 4.79E-05 3.20Eþ02 9.60E-04
1356
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
on the total sensitivity is in the following order: Injection volume, Hydraulic permeability of aquifer, Height of aquifer, Injection temperature, Anisotropy, Thermal conductivity of aquifer, Hydraulic permeability of aquitard, and Specific heat of aquifer. There is some difference in selecting important input variables between the results of this study and those obtained by Schout et al. [15]. Even though they found that variations in permeability, injection temperature and injection volume have the greatest impact, there is no change in recovery efficiency due to injection temperature in this study. Relative variations in recovery efficiency are listed in Table 9. This difference is caused by the uniformity of input parameter sampling and interaction effects. There was no consideration of interaction effects and once-at-time method was used in fractional factorial design without normalization of each factor range (without using equal sampling step). For example, thermal conductivity used for sensitivity analysis is 9 times higher than that of base case, while specific heat is 1.5 times than that of base case. In Case 2, the strongest main effects are Hydraulic permeability, Specific heat of aquifer, and Specific heat of aquitard. The sum of their main effects and first-order interaction is about 86% of the total variation. The other global sensitivity indices are extremely low and the sum of all main effects and first-order interaction is 0.8715. It means that higher order (more than second) interaction effects cover 13% of the total variation. The ranking of importance for input variables based on the total sensitivity are in the following order: Hydraulic permeability of aquifer, Specific heat of aquifer, Specific heat of aquitard, Injection temperature, Anisotropy of aquifer, Anisotropy of aquitard, Thermal conductivity of aquifer, Thermal conductivity of aquitard, Hydraulic permeability of aquitard, and Injection volume.
The gi(xi) represents the main effects and gij(xi,xj) is regarded as first-order interaction effects, and so on. Sobol [33,34] suggested the method to estimate the influence of individual variables or groups of variables on the model output. The total variance and partial variances are defined as
Z D¼
Z g 2 ðxÞdx g02
and
Di1 …is ¼
gi21 …is dxi1 …dxis
(14)
The global sensitivity indices are also defined as
Si1 …is ¼
Di1 …is D
(15)
All Si1 …is s are non-negative and their sum is
Xn Xn s¼1
S i1 < … < is i1 …is
¼1
(16)
The global sensitivity indices are usually used for ranking the variables, fixing the unessential variables, and deleting the high order members. Homma and Saltelli [35] proposed the total sensitivity index which is the total effect of each variable including its interactions with other variables.
STi ¼ Si þ Sði;iÞ ¼ 1 Si
(17)
where Si is the sum of all Si1 …is terms which do not include the i th variable. Sobol [34] also proposed the computation algorithms of Di1 …is using Monte Carlo methods. Global sensitivity indices and total sensitivity indices for Case 1 and Case 2 are listed in Table 8. The main effects of HT-ATES for Case 1 and Case 2 are plotted in Fig. 7. In Case 1, the strongest main effects are Injection volume and Hydraulic permeability of aquifer, and the sum of their main effects and first-order interaction is about 55% of the total variation. The other global sensitivity indices are extremely low and the sum of all main effects and first-order interaction is 0.549177. It means that higher order (more than second) interaction effects cover 45% of total variation. The ranking of importance for input variables based
5. Conclusion and discussion A sensitivity analysis plays an important role in presenting key parameters of a computer model dealing with the complex physical phenomena. However the reliability of sensitivity analysis depends on the number of level-combinations. In case of a computational
Table 8 Global sensitivity indices and total sensitivity indices. Case 1 S
ka
Anisotropic
Ti
Vi
H
la
kc
Ca
ka Anisotropic Ti Vi H kc Ca
1.129E-03 6.416E-12 3.178E-11 5.263E-01 2.394E-10 5.426E-13 1.116E-13 0.000Eþ00
6.416E-12 9.059E-13 6.943E-20 2.315E-11 1.197E-19 9.912E-22 1.338E-22 0.000Eþ00
3.178E-11 6.943E-20 1.709E-11 1.201E-09 1.330E-18 7.032E-21 9.392E-22 0.000Eþ00
5.263E-01 2.315E-11 1.201E-09 2.172E-02 3.430E-09 1.113E-11 6.595E-13 0.000Eþ00
2.394E-10 1.197E-19 1.330E-18 3.430E-09 2.372E-11 9.517E-21 2.733E-21 0.000Eþ00
5.426E-13 9.912E-22 7.032E-21 1.113E-11 9.517E-21 4.271E-14 5.253E-24 0.000Eþ00
1.116E-13 1.338E-22 9.392E-22 6.595E-13 2.733E-21 5.253E-24 0.000Eþ00 0.000Eþ00
0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00
ST
0.5274
3.047E-11
1.250E-09
0.5480
3.693E-09
1.171E-11
7.710E-13
0
la
Case 2 S
ka
Anisotropy
kc
Anisotropy
la
lc
Ca
Cc
Vi
Ti
ka Anisotropy kc anisotropy
Ca Cc Vi Ti
2.67E-01 5.76E-04 1.05E-10 7.72E-04 7.87E-05 2.77E-05 2.25E-01 1.23E-01 6.61E-14 3.97E-03
5.76E-04 9.12E-05 1.38E-14 6.72E-07 6.26E-08 1.73E-08 8.77E-05 1.72E-04 1.36E-16 1.29E-06
1.05E-10 1.38E-14 1.87E-12 3.48E-14 1.12E-14 5.12E-15 1.71E-11 9.58E-12 1.09E-24 1.79E-12
7.72E-04 6.72E-07 3.48E-14 1.96E-04 2.27E-07 1.51E-08 7.32E-05 1.72E-04 7.30E-17 8.30E-06
7.87E-05 6.26E-08 1.12E-14 2.27E-07 9.10E-05 2.57E-08 3.04E-05 2.50E-05 1.46E-17 6.40E-07
2.77E-05 1.73E-08 5.12E-15 1.51E-08 2.57E-08 9.81E-06 1.02E-05 1.46E-05 1.24E-17 3.39E-07
2.25E-01 8.77E-05 1.71E-11 7.32E-05 3.04E-05 1.02E-05 1.10E-01 1.05E-01 1.16E-14 1.71E-03
1.23E-01 1.72E-04 9.58E-12 1.72E-04 2.50E-05 1.46E-05 1.05E-01 2.72E-02 1.97E-14 2.78E-03
6.61E-14 1.36E-16 1.09E-24 7.30E-17 1.46E-17 1.24E-17 1.16E-14 1.97E-14 0.0Eþ00 3.18E-16
3.97E-03 1.29E-06 1.79E-12 8.30E-06 6.40E-07 3.39E-07 1.71E-03 2.78E-03 3.18E-16 3.67E-03
ST
0.6207
9.30E-04
1.35E-10
0.0012
2.26E-04
6.27E-05
0.4416
0.2577
9.80E-14
0.0121
la lc
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1357
Fig. 7. Main effect of recovery efficiency of HT-ATES.
model with many input variables, nevertheless, estimations of output variable with all level-combinations are impractical due to computation efforts. Through the design and modeling of a computer experiment, a reliable sensitivity analysis with data representing whole experiment domain is possible. In this study, the sensitivity analysis of high-temperature aquifer thermal energy
storage with single well was performed to investigate the main and interaction effects on the recovery efficiency. The design of computer experiment is conducted using the optimal Latin hypercube sampling with enhanced stochastic evolutionary based on various uniformity measurements. Then, a surrogate model, a so-called ‘metamodel’, was obtained using the Gaussian Kriging method
1358
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
ka
Ti
Vi
H
la
kc
Ca
(2013R1A2A2A01067898) and a grant (14-RDRP-B076564-01) from Regional Development Research Program funded by Ministry of Land, Infrastructure and Transport of Korean government.
49.04% 17.65%
0% 29.41%
50.96% 32.35%
0% 2.94%
0% 8.82%
0% 4.41%
0% 4.41%
References
Table 9 Relative variationa in recovery efficiency.
This study Schout et al a
Although used parameters for sensitivity analysis are different, the same parameters can be selected for comparison because the longitudinal dispersivity and anisotropy effect are negligible. In order to compare both study results, relative variation, which is defined as ratio of each variation to the sum of each variation, is adopted. Variations are computed by absolute mean value of variation in case of Schout et al., and those are obtained from global sensitivity indices in this study.
with SCAD (Smoothly Clopped Absolute Deviation) penalty. Finally, a sensitivity analysis was conducted based on a variation decomposition. This approach for sensitivity analysis for hightemperature aquifer thermal energy storage with single well allows us to compute the significance of the model parameters with representativeness, thus giving us the possibility to select the best site for HT-ATES application and realize optimum operation. The most important input variables for the Case 1 HT-ATES are injection volume with total sensitivity index of 0.5480 and hydraulic permeability of aquifer with total sensitivity index of 0.5274, and the ranking of input variables is as follows: volume, hydraulic permeability of aquifer, height of aquifer, injection temperature, anisotropy, thermal conductivity of aquifer, hydraulic permeability of aquitard, and specific heat of aquifer. Except the operation parameters such as injection volume and injection temperature, the hydraulic and thermal properties of the aquifer could be controlled by engineers. Therefore, the operation parameters are considered as key parameters for increasing the recovery efficiency. In view of this, the injection volume should be considered a priority than the injection temperature. In addition, the sensitivity analysis used in Schout et al. [15] causes misestimation of factors due to sampling and analysis problem (e.g. Relative variation of Temperature is 0% in this study, and 29.41% in study of Schout et al. [15]). Sensitivity analysis based on proper sampling method is strongly recommend. In Case 2, the most important input variables of HT-ATES are hydraulic permeability with total sensitivity index of 0.6207, specific heat of aquifer with total sensitivity index of 0.4416, and specific heat of aquitard with total sensitivity index of 0.2577 and the ranking of input variables is as follows: hydraulic permeability of aquifer, specific heat of aquifer, specific heat of aquitard, injection temperature, anisotropy of aquifer, anisotropy of aquitard, thermal conductivity of aquifer, thermal conductivity of aquitard, hydraulic permeability of aquitard, and injection volume. In controversy to Case 1, the injection temperature is more important than the injection volume. It implies that although the same computer model for a hightemperature aquifer thermal energy storage system is used, the key parameters vary with the experiment domain of hydraulic and thermal properties and the number of input variables to consider the site. This means that site-specific sensitivity analysis is needed. Therefore, utilization of key parameters is available only after a sensitivity analysis with a computer model considering all possible input parameters. Selection of a good location based on the hydraulic permeability, specific heat of aquifer and aquitard in the Gassum formation could lead a good result. Though further sensitivity analysis for HT-ATESs with several wells, which has general conditions, should be performed.
Acknowledgment This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP)
[1] Kim JK, Lee YM, Yoon WS, Jeon JS, Koo MH, Keehm YS. Numerical modeling of aquifer thermal energy storage system. Energy 2010;35:4955e65. http:// dx.doi.org/10.1016/j.energy.2010.08.029. €m G, Nordell B. Heating and cooling with UTES in [2] Andersson O, Hellstro Sweden e current situation and potential market development, FUTURESTOCK'2003. In: 9th International Conference on Thermal Energy Storage: Warsaw, Poland. Waraszawa: PW Publishing House; September 1-4, 2003. p. 207e15. [3] Godschalk MS, Bakema G. 20,000 ATES systems in the Netherlands in 2020 e major step towards a sustainable energy supply. 11/11/2009. www.iftec.es. [4] IEA. Annex VI e executive summary IF technology. The Netherlands: Arnhem; 1995. €llmann G, Hoffmann F, Bartels J. Two-year experience in the [5] Kabus F, Mo operation of an aquifer thermal energy store based on surplus heat arising from a gas and steam cogeneration plant at Neubrandenburg, NE Germany. In: Proceedings of the 10th International Conference on Thermal Energy Storage (EcoStock), New Jersey, United States; May 31 e June 2, 2006. p. 613e20. [6] Sanner B, Kabus F, Seibt P, Bartels J. Underground thermal energy storage for the German Parliament in Berlin, system concept and operational experiences. In: Proceedings World Geothermal Congress 2005, Antalya, Turkey; 24-29, April 2005. p. 8. No 1438. €i S, Driesner T. Realistic simulation of an aquifer thermal [7] Yapparova A, Mattha energy storage: effects of injection temperature, well placement and groundwater flow. Energy 2014;76:1011e8. http://dx.doi.org/10.1016/ j.energy.2014.09.018. [8] Fang KT, Li R, Sudjianto A. Design and modeling for computer experiments. 1st ed. Boca Raton, FL, USA: Chapman & Hall/CRC; 2006. [9] Jiang X, Li S, Zhang L. Sensitivity analysis of gas production from class I hydrate reservoir by depressurization. Energy 2012;39:281e5. http://dx.doi.org/ 10.1016/j.energy.2012.01.016. [10] Ashouri A, Petrini F, Bornatico R, Benz MJ. Sensitivity analysis for robust design of building energy systems. Energy 2014;76:264e75. http://dx.doi.org/ 10.1016/j.energy.2014.07.095. [11] Tarancon MA, Río PD. Assessing energy-related CO2 emissions with sensitivity analysis and input-output techniques. Energy 2012;37:161e70. http:// dx.doi.org/10.1016/j.energy.2011.07.026. [12] Miro S, Hartmann D, Schanz T. Global sensitivity analysis for subsoil parameter estimation in mechanized tunneling. Comput Geotech 2014;56:80e8. http://dx.doi.org/10.1016/j.compgeo.2013.11.003. F. Sensitivity analysis for volcanic source modeling quality assess[13] Cannavo ment and model selection. Comput Geosci 2012;44:52e9. http://dx.doi.org/ 10.1016/j.cageo.2012.03.008. [14] Rman N. Sensitivity analysis of a numerical model of sedimentary basin geothermal system. In: 30th New Zealand Geothermal Association Workshop 2008, Taupo Great Lake Centre, New Zealand; 11-12, Nov 2008. p. 298e303. [15] Schout G, Drijver B, Gutierrez-Neri M, Schotting R. Analysis of recovery efficiency in high-temperature aquifer thermal energy storage: a Rayleigh-based method. Hydrogeol J 2014;22:281e91. http://dx.doi.org/10.1007/s10040013-1050-8. [16] Bridger DW, Allen DM. Influence of geologic layering on heat transport and storage in an aquifer thermal energy storage system. Hydrogeol J 2014;22: 233e50. http://dx.doi.org/10.1007/s10040-013-1049-1. [17] Buscheck TA. The hydrothermal analysis of aquifer thermal energy storage. Berkeley, United States: University of California; 1984 [PhD Thesis]. [18] Gutierrez-Neri M, Buik N, Drikver B, Godschalk B. Analysis of recovery efficiency in a high-temperature energy storage system. In: Proceedings of the First National Congress on Geothermal Energy, Utrecht, Netherlands; 13e14, October 2011. [19] Kirkpatrick S, Gelett C, Vecchi M. Optimization by simulated annealing. Science 1983;220:671e80. [20] Jin R, Chen W, Sudjianto A. An efficient algorithm for constructing optimal design of computer experiments. J Stat Plan Inference 2005;134:268e87. http://dx.doi.org/10.1016/j.jspi.2004.02.014. [21] Fang KT, Ma CX, Winker P. Centered L2-discrepancy of random sampling and Latin hypercube design, and construction of uniform designs. Math Comp 2000;71:275e96. http://dx.doi.org/10.1090/S0025-5718-00-01281-3. [22] Hickernell FJ. A generalized discrepancy and quadrature error bound. Math Comp 1998;67:299e322. http://dx.doi.org/10.1090/S0025-5718-98-00894-1. [23] Hickernell FJ. Lattice rules: how well do they measure up? In: Hellekalek P, Larcher G, editors. Random and Quasi-random Point sets. New York: SpringerVerlag; 1998. p. 106e66. http://dx.doi.org/10.1007/978-1-4612-1702-2_3. [24] Booker AJ, Dennis JJ, Frank P, Serafini D, Torczon V, Trosset M. A rigorous framework for optimization of expensive function by surrogates. Struct Optim 1999;17:1e13. http://dx.doi.org/10.1007/BF01197708. [25] Hoffman RM, Sudjianto A, Du X, Stout J. Robust piston design and optimization using piston secondary motion analysis. Society of Automotive Engineers. 2003. http://dx.doi.org/10.4271/2003-01-0148. SAE Paper 2003-01-0148.
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359 [26] Du X, Chen W. Sequential optimization and reliability assessment method for efficient probabilistic design. ASME J Mech Des 2004;126:225e33. http:// dx.doi.org/10.1115/1.1649968. [27] Kleijnen Jack PC. Statistical tools for simulation practitioners. New York, NY, USA: Marcel Decker; 1987. [28] Simpson TW, Peplinski J, Koch PN, Allen JK. Metamodels for computer-based engineering design: survey and recommendations. Eng Comput 2001;17: 129e50. http://dx.doi.org/10.1007/PL00007198. [29] Krige DG. A statistical approach to some mine valuations and allied problems at the Witwatersrand [Master’s thesis]. University of Witwatersrand; 1951. [30] Matheron G. Principles of geostatistics. Econ Geol 1963;58:1246e66. http:// dx.doi.org/10.2113/gsecongeo.58.8.1246. [31] Li R, Sudjianto A. Analysis of computer experiments using penalized likelihood in Gaussian Kriging models. Technometrics 2005;47:111e20. http:// dx.doi.org/10.1198/004017004000000671. [32] Fan J, Li R. Variable selection via nonconcave penalized likelihood and its oracle properties. J Am Stat Assoc 2001;96:1348e60. http://dx.doi.org/ 10.1198/016214501753382273. [33] Sobol' IM. Sensitivity analysis for nonlinear mathematical models. Math Model Comput Exp 1993;1:407e14. http://dx.doi.org/10.1016/S0378-475400-00270-6. [34] Sobol' IM. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comp Simul 2001;55:271e80. http:// dx.doi.org/10.1016/S0378-4754-00-00270-6. [35] Homma T, Saltelli A. Importance measures in global sensitivity analysis of nonlinear models. Reliab Eng Syst Saf 1996;52:1e17. http://dx.doi.org/ 10.1016/0951-8320-96-00002-6.
Nomenclature ka: hydraulic permeability of aquifer (mD) kc: hydraulic permeability of aquitard (mD) H: thickness of aquifer (m) Ca: specific heat of aquifer (J/(kg K)) Cc: specific heat of aquitard (J/(kg K)) Cw: volumetric heat capacity of water (J/(m3 K)) Vp: produced volume of water (m3) Vi: injected volume of water (m3) T p : average temperature of produced water (K) Ta: ambient temperature (K) Ti: temperature of injected water (K) Greek symbols ε: recovery efficiency of HT-ATES qbi : smoothing parameter b 2 : estimate of variance s m b : estimate of mean l: regularization parameter b l: estimate of regularization parameter la: thermal conductivity of aquifer (W/(m K)) lc: Thermal conductivity of aquitard (W/(m K))
1359
Contents lists available at ScienceDirect
Energy journal homepage: www.elsevier.com/locate/energy
Sensitivity analysis of recovery efficiency in high-temperature aquifer thermal energy storage with single well Jun-Seo Jeon a, 1, Seung-Rae Lee a, *, Lisa Pasquinelli b, Ida Lykke Fabricius b a b
Department of Civil and Environmental Engineering, KAIST, Daejeon 305-701, South Korea Department of Civil Engineering, DTU BYG, Denmark
a r t i c l e i n f o
a b s t r a c t
Article history: Received 21 November 2014 Received in revised form 4 May 2015 Accepted 19 June 2015 Available online 17 July 2015
High-temperature aquifer thermal energy storage system usually shows higher performance than other borehole thermal energy storage systems. Although there is a limitation in the widespread use of the HTATES system because of several technical problems such as clogging, corrosion, etc., it is getting more attention as these issues are gradually alleviated. In this study, a sensitivity analysis of recovery efficiency in two cases of HT-ATES system with a single well is conducted to select key parameters. For a fractional factorial design used to choose input parameters with uniformity, the optimal Latin hypercube sampling with an enhanced stochastic evolutionary algorithm is considered. Then, the recovery efficiency is obtained using a computer model developed by COMSOL Multiphysics. With input and output variables, the surrogate modeling technique, namely the Gaussian-Kriging method with Smoothly Clopped Absolute Deviation Penalty, is utilized. Finally, the sensitivity analysis is performed based on the variation decomposition. According to the result of sensitivity analysis, the most important input variables are selected and confirmed to consider the interaction effects for each case and it is confirmed that key parameters vary with the experiment domain of hydraulic and thermal properties as well as the number of input variables. © 2015 Elsevier Ltd. All rights reserved.
Keywords: High-temperature aquifer thermal energy storage Sensitivity analysis Gaussian Kriging method Computational experiment Gassum formation
1. Introduction Energy storage technology is an essential technology for efficient energy management and utilization. Typical energy storage is to transform electric energy obtained by fossil fuel, nuclear power, and renewable source to mechanical, chemical energy form and then save it. It has advantages to reduce energy use, greenhouse gas emission, and cost for heating and cooling. Thermal energy storage is the storage of energy in the form of heat, and hence it may be directly utilized for industrial heat, as well as providing heating or power generation. The thermal energy storage functions to compensate for the temporal and quantitative differences in the supply and demand of energy, and energy saving through it and to improve the energy use efficiency. UTES (Underground thermal energy storage) systems are generally divided into ATES (aquifer thermal energy storage), CTES (cavern thermal energy storage), and BTES (borehole thermal
* Corresponding author. Tel.: þ82 42 350 3617; fax: þ82 42 350 7200. E-mail address: [email protected] (S.-R. Lee). 1 Current address: Department of Civil and Environmental Engineering, KAIST, 291, Gwahakro, Yuseong-gu, Daejeon 305-701, South Korea. http://dx.doi.org/10.1016/j.energy.2015.06.079 0360-5442/© 2015 Elsevier Ltd. All rights reserved.
energy storage). In ATES systems, thermal energy is stored in the ground water of an aquifer and groundwater is injected and extracted by one or more wells. CTES systems utilize underground cavern to store thermal energy using hot water and therefore thermal stratification has a key role in performance of the system. BTES systems are the most commonly used and thermal energy is stored in the storage medium, usually bedrock, using closed heat exchange pipe. The ATES system has higher system performance than the BTES system and any other systems using low temperature geothermal heat because of directly using groundwater with relatively high volumetric heat capacity [1]. Andersson et al. [2] also showed similar results through investigation on the performance factor, energy saving and payback time for ATES systems in Sweden (Table 1). In this regard, ATES systems are more suitable for large scale systems than BTES systems and mainly are used for seasonal thermal energy storage both of heating and cooling. Many applications of ATES are employed in European countries such as Netherland, Sweden and Germany. In particular, five systems in 1990's, 214 systems in 2000's, and 1300 ATES systems in 2010's so far have been registered in Netherlands. The Netherlands has a leading position in the world with rapid growth of installing ATES system [3] (Fig. 1).
1350
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
Table 1 Benefits in terms of economic and energy saving (Andersson et al. [2])*. System application
Performance factor
Energy saving (%)
Payback (year)
Direct heating and cooling HP supported heating and cooling HP supported heating only Direct cooling only
20e40
90e95
0e2
5e7
80e87
1e3
3e4
60e75
4e8
20e60
90e97
0e2
*Direct heat and cooling could be considered as HT-ATES and second and third system could be considered as LT-ATES.
ATES is generally divided into HT-ATES (high temperature ATES) system and LT-ATES (low temperature ATES) system. In case of HTATES system, the high temperature can be used directly for heating. It helps to minimize or eliminate the usage of heat pump, and then more energy is saving than LT-ATES system, which generally use heat pumps to increase the temperature. Despite the significant advantages, the HT-ATES is hard to employ in practice due to severe operating problems. It covers the clogging due to particle, gas bubbles and precipitation of minerals, and corrosion of components in the ground water system, etc [4]. According to the definition of high temperature (minimum storage loading temperatures on the order of 50 C) in ECES Annex 12, there are only two operational HT-ATES systems in Germany. These are installed in Neubrandenburg by 2004 with 80 C storage loading temperature [5] and at the Reichstag building in Berlin by 1999 with 70 C storage loading temperature [6]. However as mentioned issues are gradually alleviated, the application of HT-ATES is increasing especially in Germany. To predict and understand the behavior of sophisticated problem, numerical simulation is usually used. The simulator also can be used to assess the feasibility of the project, as well as to help in storage design and optimization of well locations, monitoring, and well steering [7]. Based on numerical simulation, sensitivity analysis is one of the main interests in decision-making processes (e.g. storage design, optimization of well locations) because it makes possible to quantify the leverage of each input variable to the output [8]. Therefore, sensitivity analysis is widely employed in the many engineering fields where dealing the complex problem
related to energy, geology, and geomechanics. For instance, sensitivity analysis is used for gas production from Class I hydrate reservoir by depressurization [9], robust design building energy systems [10], assessment of energy-related CO2 emissions [11], subsoil parameter estimation in mechanized tunneling [12], volcanic source modeling quality assessment and model selection [13], and sedimentary basin geothermal system [14]. Several studies on sensitivity analysis of ATES were also recently introduced. Schout et al. [15] suggested a method to estimate the recovery efficiency based on the Rayleigh number using results of sensitivity analysis. According to their results, performance of HT-ATES, which can be calculated in terms of recovery efficiency, depends on the hydraulic and operation properties (eg. permeability, injection temperature and injection volume). Bridger and Allen [16] investigated the influence of geologic layering on heat transport and storage in ATES. Results of the sensitivity analysis state that vertical anisotropy and hydraulic gradient are most sensitive on heat transport in a nonlayered aquifer, which is analogous to our case. Kim et al. [1] investigated the effect of injection rate, hydraulic conductivity, and distance between injection and production wells, and then they showed that system performance depends on the distance between boreholes, hydraulic conductivity of an aquifer, and the production/injection rate. Yapparova et al. [7] also performed numerical simulations to assess the effect of groundwater flow and injection temperature. Despite the results from previous several studies, the most sensitivity factor is different from ours. Even though the sensitivity analysis was used to select important factors, they did not confirm the interaction effects. Furthermore the onceat-time method was used in fractional factorial design corresponding to maximum and minimum values of each factor. It has the possibility to diminish the uniformity for sampling. In our study, we investigated main and interaction effects on the efficiency recovery of HT-ATES with single well using a proper computer experiment procedure. We considered two different cases and variations in hydraulic, thermal, and operation properties. Our findings will be helpful for the preliminary design and operation phases of the HT-ATES (high temperature-aquifer thermal energy system). This paper is structured with a flowchart illustrated in Fig. 2. First, the development of computational model for HT-ATES with single well is conducted using commercial program. Through the comparison between the recovery efficiency and tilt angle in thermal contours obtained by field experimental data and numerical simulation data, the developed model is validated.
Fig. 1. ATES systems in the Netherlands (Godschalk and Bakerma [3]).
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1351
Fig. 2. Flow chart for sensitivity analysis of HT-ATES.
Second, the fractional factorial sampling is completed based on Latin Hypercube sampling with optimization and criterion. Using that, the recovery efficiency for HT-ATES is obtained by the developed model. Third, the Gaussian Kriging model with penalized likelihood and SCAD penalty is employed to construct the metamodel. In the last, the sensitivity analysis using the Sobol's indices is performed and conclusions are presented.
2. Model construction and validation The commercial program COMSOL Multiphysics was used for modeling of the HT-ATES. The program provides a variety of physics interfaces which can be combined to handle multiphysics applications. For modeling of the HT-ATES, the fluid flow and heat transfer equations are considered as governing equations. The twodimensional axial symmetric model was considered by employing heat transfer module and Darcy's law in porous media implemented in COMSOL. The hydraulic gradient of HT-ATES is assumed to be zero since the aquifer is located below 1.5 km from the ground surface. In order to limit the number of input variables and see the change in physical properties, single injection/production well was considered. Physical and thermal properties of fluid and aquifer for the modeling are listed in Table 2. Input properties were chosen equal to those of the base case that was used in the work of Buscheck [17] for model validation later. Fluid density, viscosity, and specific heat variations resulting from the high injection water temperature into the aquifer were accounted for the simulations (Fig. 3). After studying the influence of different mesh sizes and considering the computation time, a discretized mesh with a total number of 9438 triangular elements, was adopted (Fig. 4). The main aspect controlling the feasibility of HT-ATES is the recovery efficiency. The recovery efficiency is defined as
ε¼
Vp $Cw $ T p Ta T p Ta ¼ Ti Ta Vi $Cw $ðTi Ta Þ
(1)
where V represents the amounts of water, Cw is the volumetric heat capacity of water, T is the average temperature, T is the temperature, and subscripts p, i, a represent the production, injection, and ambient. According to Gutierrez-Neri et al. [18], the recovery efficiencies were calculated based on the produced water temperature at
Table 2 Input parameter values used in the base case. Parameter
Units
Domain Aquifer
Aquitard
Horizontal permeability Vertical permeability Porosity Thermal conductivity Specific heat Ambient temperature Thickness
mD mD e W/(m K) J/(kg K) C m
15,000 1500 0.25 2.29 696 20 21
0.15 0.015 0.35 2.56 696 20 9
Parameter
Units
Operation Injection/Production
Storage and rest
Pumped volume Injection flow rate Injection temperature Duration
m3 m/s C days
55,000 5.371 105 90 90/90
e e e 90/90
fourth year operation for stabilization. Assuming 360 days in a year, 90 days of injection, 90 days of storage, 90 days of production, and 90 days of storage was considered as one cycle. Then the analysis was performed for four cycles. To validate the model, the results of Buscheck [17] and Schout et al. [15] were used. Schout et al. [15] also verified their model by comparing the angle of tilt for 14 scenarios and recovery efficiencies with those of Buscheck [17]. Mean thermal contour after first injection complete (90 days) and storage (180 days) for AL55 and AI90, which is computed by the developed model is shown in Fig. 5. The dotted line represents the angle of tilt, which is the angle of the mean temperature contour with the vertical. The average difference in the angles of tilt obtained by the developed model from those of Buscheck [17] and Schout et al. [15] are 3.45% and 2.80% (Table 3). The developed model results are in good agreement with the results found by Buscheck [17] and Schout et al. [15]. Therefore, the developed model for the ATES system using COMSOL could give a reasonable prediction of the recovery efficiencies.
3. Designing and modeling for computer experiments A general computer experiments process usually covers choosing factors and experimental domain, constructing a design of
1352
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
Fig. 3. Input parameter variation with temperature.
Fig. 4. Two-dimensional axial symmetric ATES model.
experiment matrix and metamodel searching, and interpreting the model and performing a sensitivity analysis (Fang et al. [8]). Based on the computer model, experimental domain is designated to seek which input factors contribute most to a response variable, here the recovery efficiency. In this study, two types of case were considered with different locations and input parameters. In Case 1, all input parameters and geometry are the same as those of Buscheck [17]. In Case 2, the sandstone layer of Gassum formation is considered as a potential layer for the application of ATES system in Denmark. The core specimens were used in the laboratory testing for collecting the input parameters. 8 and 10 input parameters were selected for Case 1 and Case 2, respectively, which may affect the output and the range was set from the minimum to the maximum value (Table 4). 3.1. Design of experiments domain All the level-combinations of the factors provide the most reliable design to estimate the main effects and some interaction
effects among factors. However, the number of runs for all the levelcombinations increases exponentially with the number of factors. Therefore, the fractional factorial design which is a subset of all the level-combinations is usually adopted. The selected subset should have a good representation of the whole combinations. Optimization rules listed in Table 5 have been widely used in searching for space-filling designs. Among them, the local search and threshold accepting have possibility to reach and stay in the local minimum. To address this problem, Kirkpatrick et al. [19] proposed the SA (simulated annealing) algorithm for a cost function that may possess many local minima. Jin et al. [20] also suggested the ESE (enhanced stochastic evolutionary) algorithm. Since ESE algorithm uses the two parameters which control the acceptance, it can adjust itself to suit different experimental design algorithms. SA and ESE algorithms are thus adopted to reach the global minimum in the optimization. Optimality criteria are also adopted to determine which subset is the best between the old and new designs. Fang et al. [21] proved that both star discrepancy and
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1353
Fig. 5. Mean thermal contour after first injection complete (90 days) and storage (180 days) for AL55 and AI90. With (a) kha ¼ 15D, kva ¼ 1:5D, khc ¼ 1:5D, kvc ¼ 0:15D, lc ¼ 2.56 W/ (m K), and injection temperature of 55 C, (b) kha ¼ 15D, kva ¼ 1:5D, khc ¼ 15,105 D, kvc ¼ 1:5,105 D, lc ¼ 0 W/(m K), and injection temperature of 95 C, superscripts h, v represents horizontal and vertical, subscript a, c indicates aquifer and aquitard.
Table 3 Angles of tilt obtained by Buscheck [17], Schout et al. [15] and COMSOL model. Scenario 90 Days
AI90 AT90 AL90 AI55 AT55 AL55 BI90 BT90 BL90 BI55 BT55 BL55 CI55 CT55 CL55
star L2 discrepancy are unsuitable for searching the best subset. Hickernell [22,23] suggested the centered discrepancy and wraparound discrepancy which satisfy the requirements of uniformity. Centered and wrap-around discrepancy are calculated as
180 Days
aB
aS
amodel DaBm DaSm aB
aS
amodel DaBm DaSm
46.1 40.3 36 16.6 13.9 13 71 68 66.4 44 39.2 35 71.3 67.1 67
49 43 38 13 14 14 75.5 75 73.5 47 38 33 68 66 62
48.01 1.91 45 4.7 40.31 4.31 19.02 2.42 17 3.1 14.03 1.03 75.42 4.42 73 5 68.19 1.79 47.23 3.23 43.6 4.4 48 13 73.3 2 71.35 4.25 68.19 1.19
65.5 65.5 64.5 34 34 29 84.5 85 85 66.5 66 64.5 84 84 84
66.37 70 66.25 36.52 31.7 27.75 84.28 82.6 84.2 67.7 68.19 65.77 85.71 84.57 82.87
0.99 2 2.31 6.02 3 0.03 0.08 2 5.31 0.23 5.6 15 5.3 5.35 6.19
64.6 65.5 63.9 32.2 23 24 82.8 80.9 81.7 64.5 64.9 64.1 83.7 80.8 81
1.77 4.5 2.35 4.32 8.7 3.75 1.48 1.7 2.5 3.2 3.29 1.67 2.01 3.77 1.87
0.87 4.5 1.75 2.52 2.3 1.25 0.22 2.4 0.8 1.2 2.19 1.27 1.71 0.57 1.13
ðCDðDn ÞÞ2 ¼
2 1 2 n Y s 13 2X 1 1 þ xkj 0:5 xkj 0:5 12 n 2 2 k¼1 j¼1 n X n Y s 1 X 1 1 þ 2 1 þ jxki 0:5j þ xji 0:5 2 2 n k¼1 j¼1 i¼1 1 xki xji 2 (2)
ðWDðDn ÞÞ2 ¼
s n Y s 4 1 X 3 xki xji 1 xki xji þ 2 3 2 n i¼1 k;j¼1
(3)
1354
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359 Table 6 Best subset of input variables.
Table 4 Factors and experimental domain for Case 1 and Case 2.
Case 1
Case 1 Parameters
Units
Domain
N
X1
X2
X3
X4
X5
X6
X7
X8
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.325 0.575 0.425 0.775 0.975 0.625 0.525 0.475 0.025 0.125 0.925 0.675 0.275 0.175 0.725 0.825 0.225 0.075 0.875 0.375
0.325 0.425 0.025 0.625 0.575 0.175 0.975 0.725 0.675 0.075 0.825 0.775 0.875 0.525 0.925 0.275 0.475 0.225 0.125 0.375
0.825 0.075 0.275 0.425 0.925 0.675 0.475 0.975 0.625 0.875 0.125 0.325 0.375 0.025 0.775 0.175 0.725 0.225 0.575 0.525
0.975 0.875 0.125 0.325 0.825 0.625 0.775 0.175 0.275 0.425 0.375 0.925 0.475 0.025 0.075 0.575 0.675 0.725 0.525 0.225
0.575 0.975 0.625 0.025 0.675 0.525 0.225 0.425 0.775 0.175 0.475 0.325 0.925 0.275 0.825 0.125 0.075 0.725 0.875 0.375
0.825 0.075 0.425 0.375 0.475 0.225 0.575 0.975 0.625 0.125 0.175 0.675 0.875 0.725 0.275 0.925 0.525 0.325 0.775 0.025
0.425 0.575 0.775 0.375 0.725 0.925 0.975 0.275 0.875 0.625 0.325 0.075 0.675 0.475 0.525 0.825 0.175 0.225 0.125 0.025
0.825 0.625 0.925 0.025 0.225 0.075 0.775 0.575 0.675 0.475 0.875 0.375 0.125 0.275 0.425 0.525 0.975 0.325 0.725 0.175
Aquifer
Aquitard
Horizontal permeability Anisotropy Porosity Thermal conductivity Specific heat Thickness Ambient temperature
mD e e W/(m K) J/(kg K) m C
[3000e75,000] [1e10] 0.25 [1.145e4.58] [600e1200] [10.5e42] 20
[7.5e–4~4.5e-3] 10 0.35 2.56 696 9 20
Parameters
Units
Operation Injection/Production
Storage and rest
m m/s C
[18,300e165,000] [1.65e-5e1.87e-4] [50e90]
e e e
Parameters
Units
Domain Aquifer
Aquitard
Horizontal permeability Anisotropy Porosity Thermal conductivity Specific heat Thickness Ambient temperature
mD e e W/(m K) J/(kg K) m C
[100e3000] [1e10] 0.25 [1.0e5.0] [600e1200] 21 55
[7.5e-4e4.5e-3] [1e10] 0.25 [1.0e5.0] [600e1200] 30 55
Parameters
Units
Operation
Pumped volume Injection flow rate Injection temperature
m3 m/s C
Pumped volume Injection flow rate Injection temperature
3
Case 2
Injection/Production
Storage and rest
[10,000e100,000] [1e-5e10e-5] [65e95]
e e e
Table 5 Considered optimization and criterion for Latin Hypercube sampling. Optimization
Criterion Distance
Uniformity
Threshold Accepting
Minimax/ maximin fp -criterion e e
Star discrepancy
Local Search Simulated Annealing Enhanced Stochastic Evolutionary
Star L2 discrepancy Centered discrepancy Wrap-around discrepancy
The discrepancy for the best subset selected by various optimality criteria are listed in Table 5. With any optimality, the ESE algorithm provides better uniformity than the SA algorithm. Thus the best subset is determined by the ESE with the centered discrepancy and wrap-around discrepancy. The design of experiments with n ¼ 20 and s ¼ 8 for Case 1, and n ¼ 20 and s ¼ 10 for Case 2 are listed in Table 6. The selected factors are uniformly scattered as shown in the scatter plot of LHS (20, 8) and LHS (20, 10) (Fig. 6). 3.2. Modeling of metamodel A computer model is an alternative of field and laboratory testing to simulate the complicated physical phenomena and to study any relationships among various parameters. However, the computational efforts increase dramatically for design optimization and probabilistic analysis based on the computer model. To alleviate these problems, the simpler and much faster to compute than the true model with accuracy is adopted by many researchers [24e26]. Kleijnen [27] called such an approximate model “model of
Case 2 N
X1
X2
X3
X4
X5
X6
X7
X8
X9
X10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.275 0.075 0.675 0.225 0.525 0.625 0.925 0.175 0.325 0.375 0.425 0.975 0.025 0.725 0.875 0.825 0.575 0.475 0.125 0.775
0.825 0.175 0.075 0.225 0.025 0.125 0.325 0.525 0.375 0.875 0.575 0.975 0.625 0.725 0.425 0.675 0.475 0.775 0.925 0.275
0.375 0.175 0.775 0.975 0.275 0.525 0.625 0.075 0.675 0.425 0.925 0.125 0.575 0.225 0.325 0.825 0.025 0.875 0.725 0.475
0.075 0.625 0.575 0.225 0.325 0.775 0.025 0.875 0.125 0.825 0.425 0.675 0.375 0.175 0.275 0.725 0.525 0.975 0.475 0.925
0.925 0.525 0.675 0.625 0.075 0.975 0.475 0.125 0.175 0.425 0.375 0.575 0.775 0.225 0.725 0.025 0.875 0.825 0.325 0.275
0.625 0.375 0.825 0.425 0.325 0.475 0.975 0.775 0.175 0.875 0.675 0.575 0.925 0.725 0.025 0.525 0.125 0.275 0.075 0.225
0.425 0.325 0.075 0.525 0.775 0.975 0.875 0.825 0.125 0.175 0.675 0.475 0.275 0.625 0.575 0.225 0.025 0.725 0.925 0.375
0.325 0.375 0.625 0.025 0.125 0.825 0.275 0.675 0.925 0.075 0.875 0.975 0.775 0.575 0.725 0.175 0.225 0.525 0.425 0.475
0.875 0.025 0.775 0.375 0.825 0.275 0.575 0.425 0.175 0.225 0.975 0.725 0.475 0.075 0.325 0.525 0.625 0.125 0.675 0.925
0.775 0.925 0.675 0.175 0.475 0.275 0.575 0.725 0.525 0.375 0.975 0.425 0.075 0.125 0.875 0.825 0.225 0.625 0.325 0.025
the model” or “metamodel”. Since the result of computer experiment is the determinant which does not contain any random errors, the metamodel is viewed as the regression on data with no random errors. In statistics, the regression method covers from the simple linear least squares method to the sophisticated neural network method. Simpson et al. [28] presented a survey by comparing various metamodel methods such as polynomial models, neural network method, and Kriging model. They recommended that the Kriging method is the best for a deterministic model, yet highly nonlinear with a moderate number of factors. Since most of computer models are defined based on governing equations dealing with the complex physical phenomena, they show usually nonlinear responses. Therefore, the Kriging method was selected for the construction of a metamodel in this study. After introducing the Kriging method by Krige [29], many researchers further developed his work. Matheron [30] suggested the Gaussian Kriging method for modeling spatial data in geostatistics. The Gaussian Kriging model is defined as
yðxÞ ¼
L X j¼0
bj Bj ðxÞ þ zðxÞ
(4)
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1355
Fig. 6. Scatter plot of the best subset.
where {Bj(x), j ¼ 1,…,L} is a basis over the experimental domain, and z(x) is assumed to be a Gaussian process with zero mean, variance and correlation function as given by
( rðq; s; tÞ ¼ corrðzðsÞ; zðtÞÞ ¼ exp
n X
) qk jsk tk j2
(5)
k¼1
where qk is the unknown smoothing parameter, n is the number of design variables, and sk, tk are kth kth components of the samples. With the assumption in the Kriging theory, the response function is composed of a regression model and stochastic process. Therefore, the first term of R.H.S in Eq. (4) represents the regression model which yields a global approximation of the design space and the second term of R.H.S in Eq. (4) represents the stochastic process which creates localized deviations. The predicted estimates at the untired points x can be written as
b b y ðxÞ ¼ bðxÞ b b þ r 0 ðxÞR1 ðqÞ y B b
(6) 0
where R is the correlation matrix, and r(x)¼(r(q;x1,x),,,,,r(q;xn,x)) . By maximizing the MLE (Maximum Likelihood Estimator), optimum unknown parameters, q and b could be obtained. The procedure of maximizing MLE are as follows (Fang et al. [8]): Step 1. Set the initial value of b to be (B0 B)1By, the least square estimator of b Step 2. Update (s2,q) by using
0
b b R1 ðqÞ y B b s 2 ¼ n1 y B b b
s
X Q b; s2 ; q ¼ [ b; s2 ; q þ n pl qj
where pl($) is a penalty function with a regularization parameter l. SCAD (Smoothly clopped absolute deviation) penalty was proposed by Fan and Li [32].
ðal qÞþ Iðq > lÞ for some a > 2 and q > 0 p0l ðqÞ ¼ l Iðq lÞ þ ða 1Þl (12) Even though l and a are unknown, they suggested that a ¼ 3.7 resulted in a good practical performance. Therefore, the Gaussian Kriging model based on the penalized likelihood with SCAD penalty was used in this study (Table 7). 4. Sensitivity analysis based on variation decomposition The metamodel g(x) can be represented in a functional ANOVA form.
gðxÞ ¼ g0 þ
(10)
X X gi ðxi Þ þ gij xi ; xj þ ,,, þ g1,,,s ðx1 ; …; xs Þ i
R1 0
i
gi1 ;…;it ðxi1 ; …; xit Þdxk ¼ 0, k ¼ i1,…,it.
(13)
Table 7 Parameters of Gaussian Kriging model with SCAD penalty. Parameters
(8)
(9)
(11)
j¼1
where
Step 3. Update b using
1 b B0 R1 ðqÞy b mle ¼ B0 R1 ðqÞB
Li and Sudjianto [31] proposed the penalized likelihood approach to alleviate a serious problem which occurred near the optimum value of parameter based on the MLE. The penalized likelihood is defined as
(7)
Then solve this equation by a numerical iteration algorithm such as NewtoneRaphson algorithm.
. v[ b; s2 ; q vq ¼ 0 1
n [ b; s2 ; q ¼ log s2 logjRðqÞj 2 2
0
1 b 2 y B b R1 ðqÞ y B b b 2s
Step 4. Iterate Step 2 and Step 3 until it converges
b l m b b s2 c q1 c q2 c q3 c q4 c q5 c q6 c q7 c q8 c q9 d q10
SCAD penalty Case 1
Case 2
0.44 0.614 0.026 8.42E-07 4.81E-06 5.65E-07 3.87E-04 2.93E-06 1.32E-05 3.07E-05 2.67E-23 e e
0.295 0.7663 0.0181 9.92E-06 3.20E-03 7.671 3.20E-03 7.20E-03 7.20E-03 4.79E-05 4.79E-05 3.20Eþ02 9.60E-04
1356
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
on the total sensitivity is in the following order: Injection volume, Hydraulic permeability of aquifer, Height of aquifer, Injection temperature, Anisotropy, Thermal conductivity of aquifer, Hydraulic permeability of aquitard, and Specific heat of aquifer. There is some difference in selecting important input variables between the results of this study and those obtained by Schout et al. [15]. Even though they found that variations in permeability, injection temperature and injection volume have the greatest impact, there is no change in recovery efficiency due to injection temperature in this study. Relative variations in recovery efficiency are listed in Table 9. This difference is caused by the uniformity of input parameter sampling and interaction effects. There was no consideration of interaction effects and once-at-time method was used in fractional factorial design without normalization of each factor range (without using equal sampling step). For example, thermal conductivity used for sensitivity analysis is 9 times higher than that of base case, while specific heat is 1.5 times than that of base case. In Case 2, the strongest main effects are Hydraulic permeability, Specific heat of aquifer, and Specific heat of aquitard. The sum of their main effects and first-order interaction is about 86% of the total variation. The other global sensitivity indices are extremely low and the sum of all main effects and first-order interaction is 0.8715. It means that higher order (more than second) interaction effects cover 13% of the total variation. The ranking of importance for input variables based on the total sensitivity are in the following order: Hydraulic permeability of aquifer, Specific heat of aquifer, Specific heat of aquitard, Injection temperature, Anisotropy of aquifer, Anisotropy of aquitard, Thermal conductivity of aquifer, Thermal conductivity of aquitard, Hydraulic permeability of aquitard, and Injection volume.
The gi(xi) represents the main effects and gij(xi,xj) is regarded as first-order interaction effects, and so on. Sobol [33,34] suggested the method to estimate the influence of individual variables or groups of variables on the model output. The total variance and partial variances are defined as
Z D¼
Z g 2 ðxÞdx g02
and
Di1 …is ¼
gi21 …is dxi1 …dxis
(14)
The global sensitivity indices are also defined as
Si1 …is ¼
Di1 …is D
(15)
All Si1 …is s are non-negative and their sum is
Xn Xn s¼1
S i1 < … < is i1 …is
¼1
(16)
The global sensitivity indices are usually used for ranking the variables, fixing the unessential variables, and deleting the high order members. Homma and Saltelli [35] proposed the total sensitivity index which is the total effect of each variable including its interactions with other variables.
STi ¼ Si þ Sði;iÞ ¼ 1 Si
(17)
where Si is the sum of all Si1 …is terms which do not include the i th variable. Sobol [34] also proposed the computation algorithms of Di1 …is using Monte Carlo methods. Global sensitivity indices and total sensitivity indices for Case 1 and Case 2 are listed in Table 8. The main effects of HT-ATES for Case 1 and Case 2 are plotted in Fig. 7. In Case 1, the strongest main effects are Injection volume and Hydraulic permeability of aquifer, and the sum of their main effects and first-order interaction is about 55% of the total variation. The other global sensitivity indices are extremely low and the sum of all main effects and first-order interaction is 0.549177. It means that higher order (more than second) interaction effects cover 45% of total variation. The ranking of importance for input variables based
5. Conclusion and discussion A sensitivity analysis plays an important role in presenting key parameters of a computer model dealing with the complex physical phenomena. However the reliability of sensitivity analysis depends on the number of level-combinations. In case of a computational
Table 8 Global sensitivity indices and total sensitivity indices. Case 1 S
ka
Anisotropic
Ti
Vi
H
la
kc
Ca
ka Anisotropic Ti Vi H kc Ca
1.129E-03 6.416E-12 3.178E-11 5.263E-01 2.394E-10 5.426E-13 1.116E-13 0.000Eþ00
6.416E-12 9.059E-13 6.943E-20 2.315E-11 1.197E-19 9.912E-22 1.338E-22 0.000Eþ00
3.178E-11 6.943E-20 1.709E-11 1.201E-09 1.330E-18 7.032E-21 9.392E-22 0.000Eþ00
5.263E-01 2.315E-11 1.201E-09 2.172E-02 3.430E-09 1.113E-11 6.595E-13 0.000Eþ00
2.394E-10 1.197E-19 1.330E-18 3.430E-09 2.372E-11 9.517E-21 2.733E-21 0.000Eþ00
5.426E-13 9.912E-22 7.032E-21 1.113E-11 9.517E-21 4.271E-14 5.253E-24 0.000Eþ00
1.116E-13 1.338E-22 9.392E-22 6.595E-13 2.733E-21 5.253E-24 0.000Eþ00 0.000Eþ00
0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00 0.000Eþ00
ST
0.5274
3.047E-11
1.250E-09
0.5480
3.693E-09
1.171E-11
7.710E-13
0
la
Case 2 S
ka
Anisotropy
kc
Anisotropy
la
lc
Ca
Cc
Vi
Ti
ka Anisotropy kc anisotropy
Ca Cc Vi Ti
2.67E-01 5.76E-04 1.05E-10 7.72E-04 7.87E-05 2.77E-05 2.25E-01 1.23E-01 6.61E-14 3.97E-03
5.76E-04 9.12E-05 1.38E-14 6.72E-07 6.26E-08 1.73E-08 8.77E-05 1.72E-04 1.36E-16 1.29E-06
1.05E-10 1.38E-14 1.87E-12 3.48E-14 1.12E-14 5.12E-15 1.71E-11 9.58E-12 1.09E-24 1.79E-12
7.72E-04 6.72E-07 3.48E-14 1.96E-04 2.27E-07 1.51E-08 7.32E-05 1.72E-04 7.30E-17 8.30E-06
7.87E-05 6.26E-08 1.12E-14 2.27E-07 9.10E-05 2.57E-08 3.04E-05 2.50E-05 1.46E-17 6.40E-07
2.77E-05 1.73E-08 5.12E-15 1.51E-08 2.57E-08 9.81E-06 1.02E-05 1.46E-05 1.24E-17 3.39E-07
2.25E-01 8.77E-05 1.71E-11 7.32E-05 3.04E-05 1.02E-05 1.10E-01 1.05E-01 1.16E-14 1.71E-03
1.23E-01 1.72E-04 9.58E-12 1.72E-04 2.50E-05 1.46E-05 1.05E-01 2.72E-02 1.97E-14 2.78E-03
6.61E-14 1.36E-16 1.09E-24 7.30E-17 1.46E-17 1.24E-17 1.16E-14 1.97E-14 0.0Eþ00 3.18E-16
3.97E-03 1.29E-06 1.79E-12 8.30E-06 6.40E-07 3.39E-07 1.71E-03 2.78E-03 3.18E-16 3.67E-03
ST
0.6207
9.30E-04
1.35E-10
0.0012
2.26E-04
6.27E-05
0.4416
0.2577
9.80E-14
0.0121
la lc
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
1357
Fig. 7. Main effect of recovery efficiency of HT-ATES.
model with many input variables, nevertheless, estimations of output variable with all level-combinations are impractical due to computation efforts. Through the design and modeling of a computer experiment, a reliable sensitivity analysis with data representing whole experiment domain is possible. In this study, the sensitivity analysis of high-temperature aquifer thermal energy
storage with single well was performed to investigate the main and interaction effects on the recovery efficiency. The design of computer experiment is conducted using the optimal Latin hypercube sampling with enhanced stochastic evolutionary based on various uniformity measurements. Then, a surrogate model, a so-called ‘metamodel’, was obtained using the Gaussian Kriging method
1358
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359
ka
Ti
Vi
H
la
kc
Ca
(2013R1A2A2A01067898) and a grant (14-RDRP-B076564-01) from Regional Development Research Program funded by Ministry of Land, Infrastructure and Transport of Korean government.
49.04% 17.65%
0% 29.41%
50.96% 32.35%
0% 2.94%
0% 8.82%
0% 4.41%
0% 4.41%
References
Table 9 Relative variationa in recovery efficiency.
This study Schout et al a
Although used parameters for sensitivity analysis are different, the same parameters can be selected for comparison because the longitudinal dispersivity and anisotropy effect are negligible. In order to compare both study results, relative variation, which is defined as ratio of each variation to the sum of each variation, is adopted. Variations are computed by absolute mean value of variation in case of Schout et al., and those are obtained from global sensitivity indices in this study.
with SCAD (Smoothly Clopped Absolute Deviation) penalty. Finally, a sensitivity analysis was conducted based on a variation decomposition. This approach for sensitivity analysis for hightemperature aquifer thermal energy storage with single well allows us to compute the significance of the model parameters with representativeness, thus giving us the possibility to select the best site for HT-ATES application and realize optimum operation. The most important input variables for the Case 1 HT-ATES are injection volume with total sensitivity index of 0.5480 and hydraulic permeability of aquifer with total sensitivity index of 0.5274, and the ranking of input variables is as follows: volume, hydraulic permeability of aquifer, height of aquifer, injection temperature, anisotropy, thermal conductivity of aquifer, hydraulic permeability of aquitard, and specific heat of aquifer. Except the operation parameters such as injection volume and injection temperature, the hydraulic and thermal properties of the aquifer could be controlled by engineers. Therefore, the operation parameters are considered as key parameters for increasing the recovery efficiency. In view of this, the injection volume should be considered a priority than the injection temperature. In addition, the sensitivity analysis used in Schout et al. [15] causes misestimation of factors due to sampling and analysis problem (e.g. Relative variation of Temperature is 0% in this study, and 29.41% in study of Schout et al. [15]). Sensitivity analysis based on proper sampling method is strongly recommend. In Case 2, the most important input variables of HT-ATES are hydraulic permeability with total sensitivity index of 0.6207, specific heat of aquifer with total sensitivity index of 0.4416, and specific heat of aquitard with total sensitivity index of 0.2577 and the ranking of input variables is as follows: hydraulic permeability of aquifer, specific heat of aquifer, specific heat of aquitard, injection temperature, anisotropy of aquifer, anisotropy of aquitard, thermal conductivity of aquifer, thermal conductivity of aquitard, hydraulic permeability of aquitard, and injection volume. In controversy to Case 1, the injection temperature is more important than the injection volume. It implies that although the same computer model for a hightemperature aquifer thermal energy storage system is used, the key parameters vary with the experiment domain of hydraulic and thermal properties and the number of input variables to consider the site. This means that site-specific sensitivity analysis is needed. Therefore, utilization of key parameters is available only after a sensitivity analysis with a computer model considering all possible input parameters. Selection of a good location based on the hydraulic permeability, specific heat of aquifer and aquitard in the Gassum formation could lead a good result. Though further sensitivity analysis for HT-ATESs with several wells, which has general conditions, should be performed.
Acknowledgment This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP)
[1] Kim JK, Lee YM, Yoon WS, Jeon JS, Koo MH, Keehm YS. Numerical modeling of aquifer thermal energy storage system. Energy 2010;35:4955e65. http:// dx.doi.org/10.1016/j.energy.2010.08.029. €m G, Nordell B. Heating and cooling with UTES in [2] Andersson O, Hellstro Sweden e current situation and potential market development, FUTURESTOCK'2003. In: 9th International Conference on Thermal Energy Storage: Warsaw, Poland. Waraszawa: PW Publishing House; September 1-4, 2003. p. 207e15. [3] Godschalk MS, Bakema G. 20,000 ATES systems in the Netherlands in 2020 e major step towards a sustainable energy supply. 11/11/2009. www.iftec.es. [4] IEA. Annex VI e executive summary IF technology. The Netherlands: Arnhem; 1995. €llmann G, Hoffmann F, Bartels J. Two-year experience in the [5] Kabus F, Mo operation of an aquifer thermal energy store based on surplus heat arising from a gas and steam cogeneration plant at Neubrandenburg, NE Germany. In: Proceedings of the 10th International Conference on Thermal Energy Storage (EcoStock), New Jersey, United States; May 31 e June 2, 2006. p. 613e20. [6] Sanner B, Kabus F, Seibt P, Bartels J. Underground thermal energy storage for the German Parliament in Berlin, system concept and operational experiences. In: Proceedings World Geothermal Congress 2005, Antalya, Turkey; 24-29, April 2005. p. 8. No 1438. €i S, Driesner T. Realistic simulation of an aquifer thermal [7] Yapparova A, Mattha energy storage: effects of injection temperature, well placement and groundwater flow. Energy 2014;76:1011e8. http://dx.doi.org/10.1016/ j.energy.2014.09.018. [8] Fang KT, Li R, Sudjianto A. Design and modeling for computer experiments. 1st ed. Boca Raton, FL, USA: Chapman & Hall/CRC; 2006. [9] Jiang X, Li S, Zhang L. Sensitivity analysis of gas production from class I hydrate reservoir by depressurization. Energy 2012;39:281e5. http://dx.doi.org/ 10.1016/j.energy.2012.01.016. [10] Ashouri A, Petrini F, Bornatico R, Benz MJ. Sensitivity analysis for robust design of building energy systems. Energy 2014;76:264e75. http://dx.doi.org/ 10.1016/j.energy.2014.07.095. [11] Tarancon MA, Río PD. Assessing energy-related CO2 emissions with sensitivity analysis and input-output techniques. Energy 2012;37:161e70. http:// dx.doi.org/10.1016/j.energy.2011.07.026. [12] Miro S, Hartmann D, Schanz T. Global sensitivity analysis for subsoil parameter estimation in mechanized tunneling. Comput Geotech 2014;56:80e8. http://dx.doi.org/10.1016/j.compgeo.2013.11.003. F. Sensitivity analysis for volcanic source modeling quality assess[13] Cannavo ment and model selection. Comput Geosci 2012;44:52e9. http://dx.doi.org/ 10.1016/j.cageo.2012.03.008. [14] Rman N. Sensitivity analysis of a numerical model of sedimentary basin geothermal system. In: 30th New Zealand Geothermal Association Workshop 2008, Taupo Great Lake Centre, New Zealand; 11-12, Nov 2008. p. 298e303. [15] Schout G, Drijver B, Gutierrez-Neri M, Schotting R. Analysis of recovery efficiency in high-temperature aquifer thermal energy storage: a Rayleigh-based method. Hydrogeol J 2014;22:281e91. http://dx.doi.org/10.1007/s10040013-1050-8. [16] Bridger DW, Allen DM. Influence of geologic layering on heat transport and storage in an aquifer thermal energy storage system. Hydrogeol J 2014;22: 233e50. http://dx.doi.org/10.1007/s10040-013-1049-1. [17] Buscheck TA. The hydrothermal analysis of aquifer thermal energy storage. Berkeley, United States: University of California; 1984 [PhD Thesis]. [18] Gutierrez-Neri M, Buik N, Drikver B, Godschalk B. Analysis of recovery efficiency in a high-temperature energy storage system. In: Proceedings of the First National Congress on Geothermal Energy, Utrecht, Netherlands; 13e14, October 2011. [19] Kirkpatrick S, Gelett C, Vecchi M. Optimization by simulated annealing. Science 1983;220:671e80. [20] Jin R, Chen W, Sudjianto A. An efficient algorithm for constructing optimal design of computer experiments. J Stat Plan Inference 2005;134:268e87. http://dx.doi.org/10.1016/j.jspi.2004.02.014. [21] Fang KT, Ma CX, Winker P. Centered L2-discrepancy of random sampling and Latin hypercube design, and construction of uniform designs. Math Comp 2000;71:275e96. http://dx.doi.org/10.1090/S0025-5718-00-01281-3. [22] Hickernell FJ. A generalized discrepancy and quadrature error bound. Math Comp 1998;67:299e322. http://dx.doi.org/10.1090/S0025-5718-98-00894-1. [23] Hickernell FJ. Lattice rules: how well do they measure up? In: Hellekalek P, Larcher G, editors. Random and Quasi-random Point sets. New York: SpringerVerlag; 1998. p. 106e66. http://dx.doi.org/10.1007/978-1-4612-1702-2_3. [24] Booker AJ, Dennis JJ, Frank P, Serafini D, Torczon V, Trosset M. A rigorous framework for optimization of expensive function by surrogates. Struct Optim 1999;17:1e13. http://dx.doi.org/10.1007/BF01197708. [25] Hoffman RM, Sudjianto A, Du X, Stout J. Robust piston design and optimization using piston secondary motion analysis. Society of Automotive Engineers. 2003. http://dx.doi.org/10.4271/2003-01-0148. SAE Paper 2003-01-0148.
J.-S. Jeon et al. / Energy 90 (2015) 1349e1359 [26] Du X, Chen W. Sequential optimization and reliability assessment method for efficient probabilistic design. ASME J Mech Des 2004;126:225e33. http:// dx.doi.org/10.1115/1.1649968. [27] Kleijnen Jack PC. Statistical tools for simulation practitioners. New York, NY, USA: Marcel Decker; 1987. [28] Simpson TW, Peplinski J, Koch PN, Allen JK. Metamodels for computer-based engineering design: survey and recommendations. Eng Comput 2001;17: 129e50. http://dx.doi.org/10.1007/PL00007198. [29] Krige DG. A statistical approach to some mine valuations and allied problems at the Witwatersrand [Master’s thesis]. University of Witwatersrand; 1951. [30] Matheron G. Principles of geostatistics. Econ Geol 1963;58:1246e66. http:// dx.doi.org/10.2113/gsecongeo.58.8.1246. [31] Li R, Sudjianto A. Analysis of computer experiments using penalized likelihood in Gaussian Kriging models. Technometrics 2005;47:111e20. http:// dx.doi.org/10.1198/004017004000000671. [32] Fan J, Li R. Variable selection via nonconcave penalized likelihood and its oracle properties. J Am Stat Assoc 2001;96:1348e60. http://dx.doi.org/ 10.1198/016214501753382273. [33] Sobol' IM. Sensitivity analysis for nonlinear mathematical models. Math Model Comput Exp 1993;1:407e14. http://dx.doi.org/10.1016/S0378-475400-00270-6. [34] Sobol' IM. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comp Simul 2001;55:271e80. http:// dx.doi.org/10.1016/S0378-4754-00-00270-6. [35] Homma T, Saltelli A. Importance measures in global sensitivity analysis of nonlinear models. Reliab Eng Syst Saf 1996;52:1e17. http://dx.doi.org/ 10.1016/0951-8320-96-00002-6.
Nomenclature ka: hydraulic permeability of aquifer (mD) kc: hydraulic permeability of aquitard (mD) H: thickness of aquifer (m) Ca: specific heat of aquifer (J/(kg K)) Cc: specific heat of aquitard (J/(kg K)) Cw: volumetric heat capacity of water (J/(m3 K)) Vp: produced volume of water (m3) Vi: injected volume of water (m3) T p : average temperature of produced water (K) Ta: ambient temperature (K) Ti: temperature of injected water (K) Greek symbols ε: recovery efficiency of HT-ATES qbi : smoothing parameter b 2 : estimate of variance s m b : estimate of mean l: regularization parameter b l: estimate of regularization parameter la: thermal conductivity of aquifer (W/(m K)) lc: Thermal conductivity of aquitard (W/(m K))
1359