- Email: [email protected]
Development of suitability maps for ground-coupled heat pump systems using groundwater and heat transport models
Geothermics 36 (2007) 459–472
Development of suitability maps for ground-coupled heat pump systems using groundwater and heat transport models Hikari Fujii a,∗ , Tadasuke Inatomi b , Ryuichi Itoi a , Youhei Uchida c a
Department of Earth Resources Engineering, Faculty of Engineering, Kyushu University, Motooka 744, Nishi-ku, Fukuoka 819-0395, Japan b YBM Co. Ltd., Kishiyama 589-10 Kitahata, Karatsu 847-1211, Japan c Geological Survey of Japan, AIST Tsukuba Central 7, Tsukuba 305-8567, Japan Received 7 March 2007; accepted 4 June 2007 Available online 23 August 2007
Abstract The thermophysical properties of subsurface materials (soils, sediments and rocks) and groundwater flow strongly affect the heat exchange rates of ground heat exchangers (GHEs). These rates can be maximized and the installation costs of the ground-coupled heat pump (GCHP) systems reduced by developing suitability maps based on local geological and hydrological information. Such maps were generated for the Chikushi Plain (western Japan) using field-survey data and a numerical modeling study. First, a field-wide groundwater model was developed for the area and the results matched against measured groundwater levels and vertical temperature profiles. Single GHE models were then constructed to simulate the heat exchange performance at different locations in the plain. Finally, suitability maps for GCHP systems were prepared using the results from the single GHE models. Variations in the heat exchange rates of over 40% revealed by the map were ascribed to differences in the GHE locations, confirming how important it is to use appropriate thermophysical data when designing GCHP systems. © 2007 CNR. Published by Elsevier Ltd. All rights reserved. Keywords: Geothermal heat pump; Groundwater flow; Numerical simulation; Chikushi Plain; Japan
1. Introduction Ground-coupled heat pump (GCHP) systems are considered to be one of the most energyefficient and environmentally friendly air-conditioning systems for temperate and sub-arctic zones. ∗
Corresponding author. Tel.: +81 92 802 3343; fax: +81 92 802 3343. E-mail address: [email protected] (H. Fujii).
0375-6505/$30.00 © 2007 CNR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.geothermics.2007.06.002
460
H. Fujii et al. / Geothermics 36 (2007) 459–472
In Japan, however, use of the GCHPs is still limited because of a lack of information on the advantages offered by these systems and because of their high initial cost. In some cases, these costs may be the result of an oversized design of the ground heat exchangers (GHEs), caused by a lack of reliable estimates for heat exchange rates. If these rates could be predicted with high accuracy during the design stage of GCHP systems, then we could also determine the minimum length of GHE required, and reduce the initial costs accordingly. Studies have been made on predicting the heat exchange rates at GHEs by means of thermal response tests (TRTs) and by determining the thermophysical properties of subsurface materials. Morita et al. (1984) developed a finite-difference simulation model of coaxial GHEs, while Deerman and Kavanaugh (1991) used an analytical model of GHEs based on the cylindrical source function (Ingersoll et al., 1954). Fujii et al. (2002) proposed and validated a new method for interpreting TRTs using the Horner’s plot (Horner, 1951), while Sanner et al. (2005) developed guidelines on the operation and interpretation of TRTs. Based on the above studies, the thermal performance of GHEs can therefore be predicted with reasonable accuracy, assuming that the thermophysical properties of the subsurface materials are known. So far, however, very little has been done to predict the heat exchange rates of GCHP systems on the basis of field-scale geological and hydrological data. Fujii et al. (2005) developed an approach for estimating the heat exchange performance of a large-scale GCHP system in Akita Plain, northern Japan, based on a field-wide groundwater flow and heat transport model that had been calibrated against measured data. Here, a similar, but extended, approach will be applied to Chikushi Plain, western Japan, leading to an evaluation of the heat exchange rates of existing GCHP systems and to the development of a suitability map for utilizing the low-enthalpy geothermal resources of the area. The suitability map can be used for the design of future GCHP systems at any location covered by the map. 2. Field-wide model of groundwater flow system in the Chikushi Plain Chikushi Plain is the largest alluvial plain on the island of Kyushu, one of the four major islands of Japan. It faces Ariake Bay to the south, and is surrounded by the Sefuri Mountains to the north, the Minou Mountains to the east and the Chikushi and Kishima Mountains to the west (Fig. 1). Four major rivers, the Yabe, Chikugo, Kase and Rokkaku, flow through the plain, a flat land lying only a few meters above sea level. The Rokkaku and Chikugo rivers divide the plain into three areas called, from west to east, Shiroishi, Saga and Chikugo (Fig. 1). The present study will concentrate on the Saga area, which includes the capital of Saga Prefecture, Saga City, and a major residential area. GCHP systems are expected to be installed in this part of the Chikushi Plain in the near future. Unconsolidated Quaternary deposits occur only in the central part of the Chikushi Plain, where their thickness varies with location, reaching a maximum of about 600 m. In the upper part of these deposits, the Ariake Clay (maximum thickness: 30 m) overlies the fluvial sands and gravels (and intercalated clay layers) that host the main aquifers of the plain. The basement in the study area is a Tertiary andesite. To investigate the groundwater flow system in the Chikushi Plain, Uchida et al. (2006) conducted a field survey of ground and surface waters. Electric conductivity, pH, water quality and stable isotope hydrogen and oxygen ratios were determined in surface and ground water samples. The results of their study showed that the groundwater aquifer is being recharged mainly from the mountainous regions of the Shiroishi, Saga and Chikugo areas.
H. Fujii et al. / Geothermics 36 (2007) 459–472
461
Fig. 1. Map of Chikushi Plain, Kyushu Island, western Japan.
In order to obtain the data required to construct a field-wide groundwater model for the present study, we monitored groundwater levels and downhole temperatures in Chikushi Plain wells from July 2004 to August 2005. Water levels were measured in 40 observation wells drilled to monitor land subsidence (depth range: 20–170 m), and in 22 pumped domestic wells (depth < 50 m), as shown in Fig. 2. The temperature profiles were obtained in the observation wells only since the sensors could not be run in the pumped wells because of the submersible pumps. The isotherm map at 50 m below sea level is shown in Fig. 3. No contours could be drawn in the Chikugo area because the observation wells were too shallow. The map shows that the groundwater temperatures in the central part of the plain at that elevation are quite uniform. The slightly higher values near Ariake Bay area may reflect the effects of seawater intrusion. Note that the average annual temperature of the water in Ariake Bay is about 18.6 ◦ C, rising above 20 ◦ C during the summer months. To develop a field-wide groundwater flow model of the Chikushi Plain, we used the finite element code FEFLOW (Diersch, 2002), which can numerically simulate groundwater flow and heat transport. A plan view of the computational grid covering all of Chikushi Plain and the surrounding mountain ranges is shown in Fig. 4; a 3D view is given in Fig. 5. The horizontal dimensions of the model are 70 km (E–W) and 55 km (N–S). Ariake Bay is included in the southern part of the model so as to achieve a reasonably accurate simulation of groundwater movement across the coastline.
462
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 2. Location of observation and groundwater wells in Chikushi Plain, particularly observation wells O-1 to O-4.
Fig. 3. Distribution of temperatures (in ◦ C) in the Saga area at an elevation of 50 m below sea level.
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 4. Finite-element grids used in the field-wide model. Thick line indicates the coastline.
Fig. 5. Three-dimensional view of the field-wide model.
463
464
H. Fujii et al. / Geothermics 36 (2007) 459–472
Table 1 Assumed physical properties for the different layers in the computational grid Layer
Soil/rock type
1–2 3–10
Clay Sand & gravel Andesite
11–15
Layer thickness (m) 0–15 0–75 40
Hydraulic conductivity (m/s)
Thermal conductivity [W/(m K)]
Heat capacity [J/(m3 K)]
Porosity (%)
1.0 × 10−9 1.0 × 10−5
1.3 1.6
4.92 × 106 4.92 × 106
40 15
1.0 × 10−7
3.2
4.92 × 106
5
Layers 1 and 2 of the model represent the Ariake Clay, and Layers 3–10 the fluvial Quaternary sands and gravels. Since the Ariake Clay and the fluvial Quaternary sand-and-gravel layers are found in the plain only, the thickness of Layers 1–10 differs significantly depending on location. The thickness of Layers 1 and 2 ranges from 0 to 15 m, whereas Layers 3–10 vary from 0 to 75 m, as shown in Table 1. Layers 11–15 correspond to the andesite basement that underlies the entire modelled area. Although the actual thickness of the basement is more than 200 m, for simplicity and considering that groundwater flow occurs mainly in the Quaternary deposits, in the model we set the total thickness of the basement at 200 m (i.e., five 40-m thick layers). Adopting this model of the geology, simulations of the natural-state groundwater flow in the Chikushi Plain were conducted using FEFLOW. It was assumed that the model was closed to fluid flow at its top and bottom. Rainfall infiltration was not included in the model due to a lack of measured data. Instead, we simulated the recharge of the system by fixing groundwater levels at all the outer boundaries of the model. These levels were determined as a function of surface elevation. The calculated water table levels were matched with measured data by adjusting the hydraulic conductivity of each layer as well as the groundwater levels at the outer boundaries. The simulation period was set at 1 million years, when the computed changes in groundwater levels with time became negligible.
Fig. 6. Comparison of measured and calculated groundwater levels.
H. Fujii et al. / Geothermics 36 (2007) 459–472
465
Fig. 7. Computed groundwater velocities in Layer 3. Solid black line indicates the coastline.
The comparison between computed and measured groundwater levels is given in Fig. 6; generally there is good agreement. Some discrepancies occur in wells with water levels below that of the sea. All of these wells are in the Shiroishi area, where there has been excessive groundwater pumpage for irrigation purposes over a number of decades. Since the abnormally low water tables do not reflect the field-scale groundwater flow system, the match shown in Fig. 6 is considered to be reasonably good. Based on the results of this part of the study, we obtained the hydraulic conductivities of the different layers (Table 1), as well as the groundwater velocities in Layer 3, located at the top of the aquifer (Fig. 7). Note that these velocities are generally low, less than 0.01 m/d, in the central part of the plain. Groundwater temperatures were then computed considering the results of the groundwater flow simulation, and compared against measured temperatures. The parameters used for matching purposes were the thermal conductivities of each layer and the heat flux through the bottom of the model (i.e., the basal heat flux). Considering the type of rocks and sediments used in the model and using data from the Japan Society of Thermophysical Properties (1990), we determined the heat capacity and porosity of the model layers (Table 1). Surface temperatures were obtained using the annual average temperature at Saga City (17.0 ◦ C) and assuming a decrease in ambient temperature with elevation of 7.0 ◦ C/1000 m. Outer boundary conditions were set to be adiabatic (no heat conduction), but allowing heat transfer by groundwater convection. Fig. 8 compares the measured and computed downhole temperature profiles in wells O-1 to O-4 (locations shown in Fig. 2). There is satisfactory agreement in the four wells except above 20 m depth, where ground temperatures are affected by seasonal changes in surface temperatures; there
466
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 8. Comparison of measured and calculated downhole temperature profiles.
were also reasonable fits in other observation wells. The thermal conductivities of the different model layers were determined on the basis of these natural-state thermal simulations (see Table 1); the estimated basal heat flux was 0.062 W/m2 . The results of the natural-state simulation indicate that the field-wide model successfully represented the groundwater flow system in the Chikushi Plain. The distribution of groundwater velocities and the thermophysical properties of each layer were used in developing the suitability maps (see Section 3).
H. Fujii et al. / Geothermics 36 (2007) 459–472
467
Fig. 9. Location of the eleven single GHE model sites.
3. Development of suitability maps for GCHP systems In order to prepare suitability maps for GCHP systems, single GHE models were constructed for eleven Chikushi Plain locations (Fig. 9) that might show different thermal properties and performance. Since the plain was the target area, the locations were mainly in the Saga area. Fig. 10 shows a 3D view of a 50-m long GHE model, with five 10-m thick layers. The computational mesh is refined at the center of the model to accurately reproduce the groundwater flow and heat transfer around the GHE. At each location, the thicknesses of the clay (Ariake Clay) and sand-and-gravel layers (Table 2) were based on geological field data. Note that no clays were found at Locations 1, 7 and 10, whereas clay only was encountered close to the Bay (Location 9). The groundwater velocity at each location, i.e. the weighted average of the velocities in the clay and sand-and-gravel layers, was obtained from the field-wide groundwater model (Table 2). Locations 1 and 7, near the mountainous area, show relatively high groundwater velocities (1.2–2.0 × 10−3 m/d) due to the large hydraulic gradient and the absence of the Ariake Clay. The estimated velocity at Location 9 is extremely small (2.5 × 10−8 m/d) because of the small gradient and the thick clay layer. Closed mass flow boundaries were assumed at the top and bottom of the GHE model. Hydraulic heads were defined at the lateral boundaries in order to establish the groundwater velocities obtained by the field-wide model for the corresponding locations. A constant-temperature boundary condition was applied to the top of the model, while a constant heat influx of 0.062 W/m2 , the same value as in the field-wide model, was used as the boundary conditions at the bottom. At the (groundwater) upstream lateral boundary a constant background temperature was assumed. Other
468
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 10. Three-dimensional view of the computational mesh in the single GHE model.
lateral boundaries were considered to be adiabatic, allowing heat transfer only by groundwater convection. The average temperature in Saga City, located in the center of Saga area (Fig. 1), was 5.9◦ and 28.9 ◦ C in January and August 2006, respectively, indicating that cooling loads are much Table 2 Inferred layer thicknesses and computed groundwater velocity at each location (see text) Location no.
1 2 3 4 5 6 7 8 9 10 11
Thickness (m)
Groundwater velocity (m/d)
Clay
Sand and gravel
0 12 32 17 16 3 0 18 50 0 18
50 38 18 33 34 47 50 32 0 50 32
2.0 × 10−3 2.3 × 10−4 2.0 × 10−5 1.5 × 10−4 5.3 × 10−5 2.8 × 10−4 1.2 × 10−3 3.2 × 10−4 2.5 × 10−8 2.5 × 10−4 7.7 × 10−6
H. Fujii et al. / Geothermics 36 (2007) 459–472
469
larger than heating loads in the area. At each location, heat exchange simulations were carried out assuming two different scenarios: • Case 1: 120 days of cooling per year without heating operations. • Case 2: 120 days of cooling and 90 days of heating per year. Thus, in Case 1 our GHE simulations only considered cooling (from the beginning of June to the end of September), assuming that non-geothermal energy sources are used for heating during the winter. In Case 2 a 90-day heating period (from the beginning of December to the end of February) was added to Case 1. In cooling and heating operations, the temperatures of the heat-transfer medium in the GHE are maintained at 35.0 and 5.0 ◦ C, respectively. It is assumed that double U-tubes and silica sand grout will be used when installing the GHEs. In the numerical model, the double U-tubes are treated using a triangular finite-element grid with an outer perimeter of 0.41 m, which is equal to the sum of the perimeters of the four U-tubes (O.D. = 0.033 m); this equates the surface area of the GHE in the model with that of the actual double U-tubes. The heat transfer coefficient between the heat-transfer medium and the outer face of the GHE is set at 19.7 W/(m2 K), a value that was obtained by history matching of FEFLOW simulation results against measured thermal response test data (Fujii et al., 2005). Fig. 11 shows, for Cases 1 and 2, the predicted heat exchange rates per unit length at the end of each annual cooling period (end of September) at Locations 1–3. Since it is in a foothill area, Location 1 has the largest hydraulic gradient of the three sites and, because the amount of Ariake Clay increases towards Ariake Bay, clay thickness is largest at Location 3. Thus, groundwater velocities decrease in the direction of the Bay, as shown in Table 2. Due to the differences in groundwater velocity and layer thermal conductivity, heat exchange rates vary with location. In Case 1, the rates at Locations 2 and 3 at 10 years are 12.1 and 20.2% lower than at Location 1, respectively. On the other hand, in Case 2 they decrease by 9.9 and 17.8%, respectively. The above results clearly show that groundwater velocity and the thermophysical properties of the subsurface materials strongly affect the long-term heat exchange performance of the GHEs. In Case 1, the heat exchange rates decline steeply during the 10 years of operation since the lack of a heating period, and corresponding cooling of the ground, results in a rapid increase in
Fig. 11. Predicted heat exchange rates for (a) Case 1 and (b) Case 2 (locations shown in Fig. 9).
470
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 12. Distribution of estimated heat exchange rates (in W/m) for the Saga area, Case 1.
subsurface temperatures, which lowers the efficiency of heat transfer into the ground. The heat exchange rates at 10 years are 12.1, 14.5 and 14.8% smaller than at 1 year at Locations 1, 2 and 3, respectively (illustrated in Fig. 11a). In Case 2, on the other hand, the decline in the heat exchange rate is small; the corresponding reductions in heat exchange rates are 2.7, 2.9 and 2.9% at Locations 1–3, respectively (illustrated in Fig. 11b). The suitability maps for GCHP systems (Figs. 12 and 13) show the heat exchange rates per unit length on the basis of the single GHE model at 11 locations after a 10-year operation period for Cases 1 and 2, respectively. The rates in Case 2 are more than 10% larger than in Case 1 at all locations because in Case 2 heat balance is maintained in the subsurface by ground cooling during winter heating operations. In both cases, heat exchange rates increase towards the mountains and decrease towards the coast, indicating that peripheral areas of the plain are more suitable for installing GCHP systems than the areas close to Ariake Bay. This phenomenon can be explained by the higher groundwater flow velocities and more permeable subsurface materials near the mountainous areas. The differences in heat exchange rate between Location 1 (near the mountains) and Location 9 (near the Bay) are 41.3 and 36.4% in Cases 1 and 2, respectively. To achieve the same thermal performance, a reduction of 41.3% in heat exchange rates will result in a 70% increase in the length of the GHE. Since drilling costs are one of the most important factors involved in determining the economic feasibility of a GCHP system, the availability of suitability maps should help optimize the selection of sites for future GCHP systems in a given targeted area.
H. Fujii et al. / Geothermics 36 (2007) 459–472
471
Fig. 13. Distribution of estimated heat exchange rates (in W/m) for the Saga area, Case 2.
4. Conclusions Suitability maps for the installation of ground-coupled heat pump (GCHP) systems were developed for the Chikushi Plain, western Japan, through a combination of groundwater field-survey data and numerical modeling of the groundwater flow system. A field-wide groundwater model was developed and the computed results matched with measured groundwater levels and subsurface temperature profiles. Single ground heat exchanger (GHE) models were then constructed to simulate the thermal performance of these exchangers at several locations in the plain. Suitability maps of GCHP systems were then developed by contouring the heat exchange rates obtained from the single GHE model. The suitability map showed that heat exchange rates could vary more than 40% depending on the location of the GHE because of differences in groundwater velocities and in the thermal conductivity of the subsurface formations, thus highlighting the importance of this type of map for the appropriate siting and design of GCHP systems. Acknowledgements The authors thank the Ministry of Land, Infrastructure and Transport of Japan and the Prefecture of Saga for the assistance provided during their field survey and for permission to publish the field data. This work was partly funded by Grant-in-Aid for Scientific Research (A) (16206089) and (B) (19360407) of the Japan Society for the Promotion of Science.
472
H. Fujii et al. / Geothermics 36 (2007) 459–472
References Deerman, J.D., Kavanaugh, S.P., 1991. Simulation of vertical U-tube ground-coupled heat pump systems using the cylindrical heat source solution. ASHRAE Trans. 97 (1), 287–295. Diersch, H.J.G., 2002. FEFLOW Reference Manual. Institute for Water Resources Planning and Systems Research Ltd., Berlin, Germany, 278 pp. Fujii, H., Akibayashi, S., Ohshima, K., 2002. Interpretation of thermal response tests in shallow deposits. Geothermal Res. Council Trans. 26, 143–148. Fujii, H., Itoi, R., Fujii, J., Uchida, Y., 2005. Optimizing the design of large-scale ground-coupled heat pump systems using groundwater and heat transport modeling. Geothermics 34, 347–364. Horner, D.R., 1951. Pressure build-up in wells. In: Proceedings of the Third World Petroleum Congress, The Hague, The Netherlands, Sec. 11, pp. 503–523. Ingersoll, L.R., Zobel, O.J., Ingersoll, A.C., 1954. Heat Conduction with Engineering, Geological, and Other Applications. McGraw-Hill, New York, NY, USA, 325 pp. Japan Society of Thermophysical Properties, 1990. Thermophysical Properties Handbook. Yokendo Co., Tokyo, Japan, 489 pp. Morita, K., Yamaguchi, T., Karasawa, H., Hayamizu, H., 1984. Development and evaluation of temperature simulation code for geothermal wells. J. Mining Metall. Inst. Jpn. 100, 1045–1051 (in Japanese with English abstract). Sanner, B., Hellstr¨om, G., Spitler, J., Gehlin, S., 2005. Thermal response test—current status and world-wide application. In: Proceedings World Geothermal Congress 2005, Paper 1436, Antalya, Turkey, 24–29 April, 9 pp. Uchida, Y., Inatomi, T., Fujii, H., 2006. Characteristics of water quality, stable isotope ratios of hydrogen and oxygen, and subsurface thermal gradient in the Chikushi Plain. J. Japanese Assoc. Hydrol. Sci. 36, 197–204 (in Japanese with English abstract).
Development of suitability maps for ground-coupled heat pump systems using groundwater and heat transport models Hikari Fujii a,∗ , Tadasuke Inatomi b , Ryuichi Itoi a , Youhei Uchida c a
Department of Earth Resources Engineering, Faculty of Engineering, Kyushu University, Motooka 744, Nishi-ku, Fukuoka 819-0395, Japan b YBM Co. Ltd., Kishiyama 589-10 Kitahata, Karatsu 847-1211, Japan c Geological Survey of Japan, AIST Tsukuba Central 7, Tsukuba 305-8567, Japan Received 7 March 2007; accepted 4 June 2007 Available online 23 August 2007
Abstract The thermophysical properties of subsurface materials (soils, sediments and rocks) and groundwater flow strongly affect the heat exchange rates of ground heat exchangers (GHEs). These rates can be maximized and the installation costs of the ground-coupled heat pump (GCHP) systems reduced by developing suitability maps based on local geological and hydrological information. Such maps were generated for the Chikushi Plain (western Japan) using field-survey data and a numerical modeling study. First, a field-wide groundwater model was developed for the area and the results matched against measured groundwater levels and vertical temperature profiles. Single GHE models were then constructed to simulate the heat exchange performance at different locations in the plain. Finally, suitability maps for GCHP systems were prepared using the results from the single GHE models. Variations in the heat exchange rates of over 40% revealed by the map were ascribed to differences in the GHE locations, confirming how important it is to use appropriate thermophysical data when designing GCHP systems. © 2007 CNR. Published by Elsevier Ltd. All rights reserved. Keywords: Geothermal heat pump; Groundwater flow; Numerical simulation; Chikushi Plain; Japan
1. Introduction Ground-coupled heat pump (GCHP) systems are considered to be one of the most energyefficient and environmentally friendly air-conditioning systems for temperate and sub-arctic zones. ∗
Corresponding author. Tel.: +81 92 802 3343; fax: +81 92 802 3343. E-mail address: [email protected] (H. Fujii).
0375-6505/$30.00 © 2007 CNR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.geothermics.2007.06.002
460
H. Fujii et al. / Geothermics 36 (2007) 459–472
In Japan, however, use of the GCHPs is still limited because of a lack of information on the advantages offered by these systems and because of their high initial cost. In some cases, these costs may be the result of an oversized design of the ground heat exchangers (GHEs), caused by a lack of reliable estimates for heat exchange rates. If these rates could be predicted with high accuracy during the design stage of GCHP systems, then we could also determine the minimum length of GHE required, and reduce the initial costs accordingly. Studies have been made on predicting the heat exchange rates at GHEs by means of thermal response tests (TRTs) and by determining the thermophysical properties of subsurface materials. Morita et al. (1984) developed a finite-difference simulation model of coaxial GHEs, while Deerman and Kavanaugh (1991) used an analytical model of GHEs based on the cylindrical source function (Ingersoll et al., 1954). Fujii et al. (2002) proposed and validated a new method for interpreting TRTs using the Horner’s plot (Horner, 1951), while Sanner et al. (2005) developed guidelines on the operation and interpretation of TRTs. Based on the above studies, the thermal performance of GHEs can therefore be predicted with reasonable accuracy, assuming that the thermophysical properties of the subsurface materials are known. So far, however, very little has been done to predict the heat exchange rates of GCHP systems on the basis of field-scale geological and hydrological data. Fujii et al. (2005) developed an approach for estimating the heat exchange performance of a large-scale GCHP system in Akita Plain, northern Japan, based on a field-wide groundwater flow and heat transport model that had been calibrated against measured data. Here, a similar, but extended, approach will be applied to Chikushi Plain, western Japan, leading to an evaluation of the heat exchange rates of existing GCHP systems and to the development of a suitability map for utilizing the low-enthalpy geothermal resources of the area. The suitability map can be used for the design of future GCHP systems at any location covered by the map. 2. Field-wide model of groundwater flow system in the Chikushi Plain Chikushi Plain is the largest alluvial plain on the island of Kyushu, one of the four major islands of Japan. It faces Ariake Bay to the south, and is surrounded by the Sefuri Mountains to the north, the Minou Mountains to the east and the Chikushi and Kishima Mountains to the west (Fig. 1). Four major rivers, the Yabe, Chikugo, Kase and Rokkaku, flow through the plain, a flat land lying only a few meters above sea level. The Rokkaku and Chikugo rivers divide the plain into three areas called, from west to east, Shiroishi, Saga and Chikugo (Fig. 1). The present study will concentrate on the Saga area, which includes the capital of Saga Prefecture, Saga City, and a major residential area. GCHP systems are expected to be installed in this part of the Chikushi Plain in the near future. Unconsolidated Quaternary deposits occur only in the central part of the Chikushi Plain, where their thickness varies with location, reaching a maximum of about 600 m. In the upper part of these deposits, the Ariake Clay (maximum thickness: 30 m) overlies the fluvial sands and gravels (and intercalated clay layers) that host the main aquifers of the plain. The basement in the study area is a Tertiary andesite. To investigate the groundwater flow system in the Chikushi Plain, Uchida et al. (2006) conducted a field survey of ground and surface waters. Electric conductivity, pH, water quality and stable isotope hydrogen and oxygen ratios were determined in surface and ground water samples. The results of their study showed that the groundwater aquifer is being recharged mainly from the mountainous regions of the Shiroishi, Saga and Chikugo areas.
H. Fujii et al. / Geothermics 36 (2007) 459–472
461
Fig. 1. Map of Chikushi Plain, Kyushu Island, western Japan.
In order to obtain the data required to construct a field-wide groundwater model for the present study, we monitored groundwater levels and downhole temperatures in Chikushi Plain wells from July 2004 to August 2005. Water levels were measured in 40 observation wells drilled to monitor land subsidence (depth range: 20–170 m), and in 22 pumped domestic wells (depth < 50 m), as shown in Fig. 2. The temperature profiles were obtained in the observation wells only since the sensors could not be run in the pumped wells because of the submersible pumps. The isotherm map at 50 m below sea level is shown in Fig. 3. No contours could be drawn in the Chikugo area because the observation wells were too shallow. The map shows that the groundwater temperatures in the central part of the plain at that elevation are quite uniform. The slightly higher values near Ariake Bay area may reflect the effects of seawater intrusion. Note that the average annual temperature of the water in Ariake Bay is about 18.6 ◦ C, rising above 20 ◦ C during the summer months. To develop a field-wide groundwater flow model of the Chikushi Plain, we used the finite element code FEFLOW (Diersch, 2002), which can numerically simulate groundwater flow and heat transport. A plan view of the computational grid covering all of Chikushi Plain and the surrounding mountain ranges is shown in Fig. 4; a 3D view is given in Fig. 5. The horizontal dimensions of the model are 70 km (E–W) and 55 km (N–S). Ariake Bay is included in the southern part of the model so as to achieve a reasonably accurate simulation of groundwater movement across the coastline.
462
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 2. Location of observation and groundwater wells in Chikushi Plain, particularly observation wells O-1 to O-4.
Fig. 3. Distribution of temperatures (in ◦ C) in the Saga area at an elevation of 50 m below sea level.
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 4. Finite-element grids used in the field-wide model. Thick line indicates the coastline.
Fig. 5. Three-dimensional view of the field-wide model.
463
464
H. Fujii et al. / Geothermics 36 (2007) 459–472
Table 1 Assumed physical properties for the different layers in the computational grid Layer
Soil/rock type
1–2 3–10
Clay Sand & gravel Andesite
11–15
Layer thickness (m) 0–15 0–75 40
Hydraulic conductivity (m/s)
Thermal conductivity [W/(m K)]
Heat capacity [J/(m3 K)]
Porosity (%)
1.0 × 10−9 1.0 × 10−5
1.3 1.6
4.92 × 106 4.92 × 106
40 15
1.0 × 10−7
3.2
4.92 × 106
5
Layers 1 and 2 of the model represent the Ariake Clay, and Layers 3–10 the fluvial Quaternary sands and gravels. Since the Ariake Clay and the fluvial Quaternary sand-and-gravel layers are found in the plain only, the thickness of Layers 1–10 differs significantly depending on location. The thickness of Layers 1 and 2 ranges from 0 to 15 m, whereas Layers 3–10 vary from 0 to 75 m, as shown in Table 1. Layers 11–15 correspond to the andesite basement that underlies the entire modelled area. Although the actual thickness of the basement is more than 200 m, for simplicity and considering that groundwater flow occurs mainly in the Quaternary deposits, in the model we set the total thickness of the basement at 200 m (i.e., five 40-m thick layers). Adopting this model of the geology, simulations of the natural-state groundwater flow in the Chikushi Plain were conducted using FEFLOW. It was assumed that the model was closed to fluid flow at its top and bottom. Rainfall infiltration was not included in the model due to a lack of measured data. Instead, we simulated the recharge of the system by fixing groundwater levels at all the outer boundaries of the model. These levels were determined as a function of surface elevation. The calculated water table levels were matched with measured data by adjusting the hydraulic conductivity of each layer as well as the groundwater levels at the outer boundaries. The simulation period was set at 1 million years, when the computed changes in groundwater levels with time became negligible.
Fig. 6. Comparison of measured and calculated groundwater levels.
H. Fujii et al. / Geothermics 36 (2007) 459–472
465
Fig. 7. Computed groundwater velocities in Layer 3. Solid black line indicates the coastline.
The comparison between computed and measured groundwater levels is given in Fig. 6; generally there is good agreement. Some discrepancies occur in wells with water levels below that of the sea. All of these wells are in the Shiroishi area, where there has been excessive groundwater pumpage for irrigation purposes over a number of decades. Since the abnormally low water tables do not reflect the field-scale groundwater flow system, the match shown in Fig. 6 is considered to be reasonably good. Based on the results of this part of the study, we obtained the hydraulic conductivities of the different layers (Table 1), as well as the groundwater velocities in Layer 3, located at the top of the aquifer (Fig. 7). Note that these velocities are generally low, less than 0.01 m/d, in the central part of the plain. Groundwater temperatures were then computed considering the results of the groundwater flow simulation, and compared against measured temperatures. The parameters used for matching purposes were the thermal conductivities of each layer and the heat flux through the bottom of the model (i.e., the basal heat flux). Considering the type of rocks and sediments used in the model and using data from the Japan Society of Thermophysical Properties (1990), we determined the heat capacity and porosity of the model layers (Table 1). Surface temperatures were obtained using the annual average temperature at Saga City (17.0 ◦ C) and assuming a decrease in ambient temperature with elevation of 7.0 ◦ C/1000 m. Outer boundary conditions were set to be adiabatic (no heat conduction), but allowing heat transfer by groundwater convection. Fig. 8 compares the measured and computed downhole temperature profiles in wells O-1 to O-4 (locations shown in Fig. 2). There is satisfactory agreement in the four wells except above 20 m depth, where ground temperatures are affected by seasonal changes in surface temperatures; there
466
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 8. Comparison of measured and calculated downhole temperature profiles.
were also reasonable fits in other observation wells. The thermal conductivities of the different model layers were determined on the basis of these natural-state thermal simulations (see Table 1); the estimated basal heat flux was 0.062 W/m2 . The results of the natural-state simulation indicate that the field-wide model successfully represented the groundwater flow system in the Chikushi Plain. The distribution of groundwater velocities and the thermophysical properties of each layer were used in developing the suitability maps (see Section 3).
H. Fujii et al. / Geothermics 36 (2007) 459–472
467
Fig. 9. Location of the eleven single GHE model sites.
3. Development of suitability maps for GCHP systems In order to prepare suitability maps for GCHP systems, single GHE models were constructed for eleven Chikushi Plain locations (Fig. 9) that might show different thermal properties and performance. Since the plain was the target area, the locations were mainly in the Saga area. Fig. 10 shows a 3D view of a 50-m long GHE model, with five 10-m thick layers. The computational mesh is refined at the center of the model to accurately reproduce the groundwater flow and heat transfer around the GHE. At each location, the thicknesses of the clay (Ariake Clay) and sand-and-gravel layers (Table 2) were based on geological field data. Note that no clays were found at Locations 1, 7 and 10, whereas clay only was encountered close to the Bay (Location 9). The groundwater velocity at each location, i.e. the weighted average of the velocities in the clay and sand-and-gravel layers, was obtained from the field-wide groundwater model (Table 2). Locations 1 and 7, near the mountainous area, show relatively high groundwater velocities (1.2–2.0 × 10−3 m/d) due to the large hydraulic gradient and the absence of the Ariake Clay. The estimated velocity at Location 9 is extremely small (2.5 × 10−8 m/d) because of the small gradient and the thick clay layer. Closed mass flow boundaries were assumed at the top and bottom of the GHE model. Hydraulic heads were defined at the lateral boundaries in order to establish the groundwater velocities obtained by the field-wide model for the corresponding locations. A constant-temperature boundary condition was applied to the top of the model, while a constant heat influx of 0.062 W/m2 , the same value as in the field-wide model, was used as the boundary conditions at the bottom. At the (groundwater) upstream lateral boundary a constant background temperature was assumed. Other
468
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 10. Three-dimensional view of the computational mesh in the single GHE model.
lateral boundaries were considered to be adiabatic, allowing heat transfer only by groundwater convection. The average temperature in Saga City, located in the center of Saga area (Fig. 1), was 5.9◦ and 28.9 ◦ C in January and August 2006, respectively, indicating that cooling loads are much Table 2 Inferred layer thicknesses and computed groundwater velocity at each location (see text) Location no.
1 2 3 4 5 6 7 8 9 10 11
Thickness (m)
Groundwater velocity (m/d)
Clay
Sand and gravel
0 12 32 17 16 3 0 18 50 0 18
50 38 18 33 34 47 50 32 0 50 32
2.0 × 10−3 2.3 × 10−4 2.0 × 10−5 1.5 × 10−4 5.3 × 10−5 2.8 × 10−4 1.2 × 10−3 3.2 × 10−4 2.5 × 10−8 2.5 × 10−4 7.7 × 10−6
H. Fujii et al. / Geothermics 36 (2007) 459–472
469
larger than heating loads in the area. At each location, heat exchange simulations were carried out assuming two different scenarios: • Case 1: 120 days of cooling per year without heating operations. • Case 2: 120 days of cooling and 90 days of heating per year. Thus, in Case 1 our GHE simulations only considered cooling (from the beginning of June to the end of September), assuming that non-geothermal energy sources are used for heating during the winter. In Case 2 a 90-day heating period (from the beginning of December to the end of February) was added to Case 1. In cooling and heating operations, the temperatures of the heat-transfer medium in the GHE are maintained at 35.0 and 5.0 ◦ C, respectively. It is assumed that double U-tubes and silica sand grout will be used when installing the GHEs. In the numerical model, the double U-tubes are treated using a triangular finite-element grid with an outer perimeter of 0.41 m, which is equal to the sum of the perimeters of the four U-tubes (O.D. = 0.033 m); this equates the surface area of the GHE in the model with that of the actual double U-tubes. The heat transfer coefficient between the heat-transfer medium and the outer face of the GHE is set at 19.7 W/(m2 K), a value that was obtained by history matching of FEFLOW simulation results against measured thermal response test data (Fujii et al., 2005). Fig. 11 shows, for Cases 1 and 2, the predicted heat exchange rates per unit length at the end of each annual cooling period (end of September) at Locations 1–3. Since it is in a foothill area, Location 1 has the largest hydraulic gradient of the three sites and, because the amount of Ariake Clay increases towards Ariake Bay, clay thickness is largest at Location 3. Thus, groundwater velocities decrease in the direction of the Bay, as shown in Table 2. Due to the differences in groundwater velocity and layer thermal conductivity, heat exchange rates vary with location. In Case 1, the rates at Locations 2 and 3 at 10 years are 12.1 and 20.2% lower than at Location 1, respectively. On the other hand, in Case 2 they decrease by 9.9 and 17.8%, respectively. The above results clearly show that groundwater velocity and the thermophysical properties of the subsurface materials strongly affect the long-term heat exchange performance of the GHEs. In Case 1, the heat exchange rates decline steeply during the 10 years of operation since the lack of a heating period, and corresponding cooling of the ground, results in a rapid increase in
Fig. 11. Predicted heat exchange rates for (a) Case 1 and (b) Case 2 (locations shown in Fig. 9).
470
H. Fujii et al. / Geothermics 36 (2007) 459–472
Fig. 12. Distribution of estimated heat exchange rates (in W/m) for the Saga area, Case 1.
subsurface temperatures, which lowers the efficiency of heat transfer into the ground. The heat exchange rates at 10 years are 12.1, 14.5 and 14.8% smaller than at 1 year at Locations 1, 2 and 3, respectively (illustrated in Fig. 11a). In Case 2, on the other hand, the decline in the heat exchange rate is small; the corresponding reductions in heat exchange rates are 2.7, 2.9 and 2.9% at Locations 1–3, respectively (illustrated in Fig. 11b). The suitability maps for GCHP systems (Figs. 12 and 13) show the heat exchange rates per unit length on the basis of the single GHE model at 11 locations after a 10-year operation period for Cases 1 and 2, respectively. The rates in Case 2 are more than 10% larger than in Case 1 at all locations because in Case 2 heat balance is maintained in the subsurface by ground cooling during winter heating operations. In both cases, heat exchange rates increase towards the mountains and decrease towards the coast, indicating that peripheral areas of the plain are more suitable for installing GCHP systems than the areas close to Ariake Bay. This phenomenon can be explained by the higher groundwater flow velocities and more permeable subsurface materials near the mountainous areas. The differences in heat exchange rate between Location 1 (near the mountains) and Location 9 (near the Bay) are 41.3 and 36.4% in Cases 1 and 2, respectively. To achieve the same thermal performance, a reduction of 41.3% in heat exchange rates will result in a 70% increase in the length of the GHE. Since drilling costs are one of the most important factors involved in determining the economic feasibility of a GCHP system, the availability of suitability maps should help optimize the selection of sites for future GCHP systems in a given targeted area.
H. Fujii et al. / Geothermics 36 (2007) 459–472
471
Fig. 13. Distribution of estimated heat exchange rates (in W/m) for the Saga area, Case 2.
4. Conclusions Suitability maps for the installation of ground-coupled heat pump (GCHP) systems were developed for the Chikushi Plain, western Japan, through a combination of groundwater field-survey data and numerical modeling of the groundwater flow system. A field-wide groundwater model was developed and the computed results matched with measured groundwater levels and subsurface temperature profiles. Single ground heat exchanger (GHE) models were then constructed to simulate the thermal performance of these exchangers at several locations in the plain. Suitability maps of GCHP systems were then developed by contouring the heat exchange rates obtained from the single GHE model. The suitability map showed that heat exchange rates could vary more than 40% depending on the location of the GHE because of differences in groundwater velocities and in the thermal conductivity of the subsurface formations, thus highlighting the importance of this type of map for the appropriate siting and design of GCHP systems. Acknowledgements The authors thank the Ministry of Land, Infrastructure and Transport of Japan and the Prefecture of Saga for the assistance provided during their field survey and for permission to publish the field data. This work was partly funded by Grant-in-Aid for Scientific Research (A) (16206089) and (B) (19360407) of the Japan Society for the Promotion of Science.
472
H. Fujii et al. / Geothermics 36 (2007) 459–472
References Deerman, J.D., Kavanaugh, S.P., 1991. Simulation of vertical U-tube ground-coupled heat pump systems using the cylindrical heat source solution. ASHRAE Trans. 97 (1), 287–295. Diersch, H.J.G., 2002. FEFLOW Reference Manual. Institute for Water Resources Planning and Systems Research Ltd., Berlin, Germany, 278 pp. Fujii, H., Akibayashi, S., Ohshima, K., 2002. Interpretation of thermal response tests in shallow deposits. Geothermal Res. Council Trans. 26, 143–148. Fujii, H., Itoi, R., Fujii, J., Uchida, Y., 2005. Optimizing the design of large-scale ground-coupled heat pump systems using groundwater and heat transport modeling. Geothermics 34, 347–364. Horner, D.R., 1951. Pressure build-up in wells. In: Proceedings of the Third World Petroleum Congress, The Hague, The Netherlands, Sec. 11, pp. 503–523. Ingersoll, L.R., Zobel, O.J., Ingersoll, A.C., 1954. Heat Conduction with Engineering, Geological, and Other Applications. McGraw-Hill, New York, NY, USA, 325 pp. Japan Society of Thermophysical Properties, 1990. Thermophysical Properties Handbook. Yokendo Co., Tokyo, Japan, 489 pp. Morita, K., Yamaguchi, T., Karasawa, H., Hayamizu, H., 1984. Development and evaluation of temperature simulation code for geothermal wells. J. Mining Metall. Inst. Jpn. 100, 1045–1051 (in Japanese with English abstract). Sanner, B., Hellstr¨om, G., Spitler, J., Gehlin, S., 2005. Thermal response test—current status and world-wide application. In: Proceedings World Geothermal Congress 2005, Paper 1436, Antalya, Turkey, 24–29 April, 9 pp. Uchida, Y., Inatomi, T., Fujii, H., 2006. Characteristics of water quality, stable isotope ratios of hydrogen and oxygen, and subsurface thermal gradient in the Chikushi Plain. J. Japanese Assoc. Hydrol. Sci. 36, 197–204 (in Japanese with English abstract).