Assessment of extreme hydrological conditions in the Bothnian Bay, Baltic Sea, and the impact of the nuclear power plant “Hanhikivi-1” on the local thermal regime

Abstract

The results of the study aimed to assess the influence of future nuclear power plant “Hanhikivi-1” upon the local thermal conditions in the Bothnian Bay in the Baltic Sea are presented. A number of experiments with different numerical models were also carried out in order to estimate the extreme hydro-meteorological conditions in the area of the construction. The numerical experiments were fulfilled both with analytically specified external forcing and with real external forcing for 2 years: a cold year (2010) and a warm year (2014). The study has shown that the extreme values of sea level and water temperature and the characteristics of wind waves and sea ice in the vicinity of the future nuclear power plant can be significant and sometimes catastrophic. Permanent release of heat into the marine environment from an operating nuclear power plant will lead to a strong increase in temperature and the disappearance of ice cover within a 2 km vicinity of the station. These effects should be taken into account when assessing local climate changes in the future.

Introduction

In recent decades the use of nuclear energy has been extended to a large scale. New nuclear power plants (NPPs) are designed and constructed, including those situated on the shores of seas and oceans, which provides free access to water needed for cooling processes, as discussed in Rubbelke and Vogele (2010). However, during the construction of these power plants, it is absolutely necessary to carry out a preliminary examination, including the assessment of risks associated with extreme natural conditions that may lead to technogenic disasters, as recently happened at the Japanese NPP in Fukushima, which was damaged during the earthquake and subsequent tsunami (see Acton and Hibbs, 2012; Buesseler et al., 2011; Srinivasan and Gopi Rethinaraj, 2013). There exists a twofold oppositely directed influence: (1) NPPs are affected by the environment and (2) NPPs impact the environment with possible negative effects manifested, in particular, in the release of warm cooling water (Chuang et al., 2009; Thermal standards for cooling water …, 2011).

Ensuring the safety of the operation of existing and planned NPPs requires solving the following two scientific problems: (1) evaluation of extreme external conditions (meteorological, hydrological, seismic) followed by an assessment of their impact on the NPPs and (2) carrying out an environmental impact assessment (EIA) for NPPs. In contrast to the first problem, which has been studied extensively in meteorology and oceanography in recent decades, researches on EIA methods, especially for the marine environment in the case of existing and planned NPPs operating in a regular mode, are not numerous (see, for example, Zeng et al., 2002; Abbaspour et al., 2005; Kaplan et al., 2016).

The environmental impact of a NPP results from the nuclear fuel cycle, the effects of nuclear accidents, and NPP operation. It is known that the greenhouse gas emissions from nuclear fission power are much smaller than those associated with coal, oil, and gas, and the routine health risks are much smaller than those associated with coal. However, there is a “catastrophic risk” potential if containment fails (von Hippel, 2010), which in nuclear reactors can be brought about by over-heated fuels melting and releasing large quantities of fission products into the environment. This potential risk could wipe out the benefits. Some predictions of the impact of severe NPP accidents on radionuclide contamination of the near-surface environment are given in Rumynin (2015).

The environmental impact of NPPs due to operation has been studied less than nuclear fission effects. During the process of nuclear power generation, large volumes of water are used. The uranium fuel inside reactors undergoes induced nuclear fission, which releases great amounts of energy that is used to heat water. The water turns into steam and rotates a turbine, creating electricity. Nuclear plants must collect around 600 gallons (2.27 m${}^{\mathrm{3}}$) of water per megawatt hour (MW h) for this process (Tellinghuisen, 2016); thus, the plants are built near bodies of water.

As with all thermoelectric plants, NPPs need cooling systems. The most common systems for thermal power plants, including nuclear, are

• once-through cooling, in which water is drawn from a large water body, passes through the cooling system, and then flows back into the water body;

• cooling pond, in which water is drawn from a pond dedicated to the purpose, passes through the cooling system, and then returns to the pond;

• cooling towers, in which water recirculates through the cooling system until it evaporates from the tower.

Nuclear plants exchange 60 to 70 % of their thermal energy by cycling with a body of water or by the evaporation of water through a cooling tower. According to World Nuclear Association data (http://www.world-nuclear.org), this thermal efficiency is somewhat lower than that of coal-fired power plants, thus creating more waste heat.

When taking in water for the cooling process, nuclear plants, like all thermal power plants, use special structures. Water is often drawn through screens to minimize the entry of debris. The problem is that many aquatic organisms are trapped and killed against the screens, through a process known as impingement. Aquatic organisms small enough to pass through the screens are subject to toxic stress in a process known as entrainment. Billions of marine organisms, such as fish, seals, shellfish, and turtles, which are essential to the food chain, are sucked into the cooling systems and destroyed.

The long-term experience related to the worldwide operation of NPPs shows that, under normal safe operating conditions, the non-radiological impact on the environment becomes dominant. One of the major factors is the heat pollution of the surface water bodies due to the discharge of waste heat from the condensers of NPPs. If heated condenser water is not cooled for reuse in a cooling tower, the waste heat may be discharged into either artificial reservoirs (ponds) or directly into surface waters like rivers, lakes, and sea bays.

NPPs release warm water into the sea, which can significantly affect the functioning of marine ecosystems on a local scale. In Hong and Guixiang (2012), a significant negative impact of the warm cooling water ejection of the coastal power plant located in the tidal Xiangshan Bay (China) on phytoplankton was shown. Similar conclusions are given in Chuang et al. (2009) after the assessment of the impact of discharged water from a coastal nuclear power plant located in Taiwan. Also, the impact assessment of the discharged water from a NPP located on the Atlantic coast in Brazil (Ilha Grande Bay) has shown the changes in the composition and structure of marine fish species as it was shown by Teixeira et al. (2009). All the studies were carried out on the basis of field measurements and observations.

In the work of Bork and Maier-Reimer (1978), by means of numerical simulation, the thermal regime in the tidal river Elbe was reproduced. As was expected, the results showed a clear oscillatory nature of the spread of warm water induced by the tidal currents. It can be assumed that, in contrast to the spread of heated water in a tidal river that always has its own flow velocity, in tidal coastal areas the oscillating contribution would be more pronounced. In Abbaspour et al. (2005), modeling of the warm waters spread from the coastal thermal power plant (Bandar Abbas thermal power plant, BATP) located in the Persian Gulf was carried out and shows good results of this method in the prediction of the discharged water spread in the basin with strong tidal oscillations. In the work of Zeng et al. (2002), the physical and numerical modeling of the warm discharged water transport spreading from a coastal nuclear power plant located in the tidal area of Daya Bay (China) near Hong Kong was carried out and also performed well in the simulation of the studied process.

In addition to field observations and numerical modeling, the usage of satellite data can also be useful in such assessments, as it was shown, for example, in the work of Chen et al. (2003), which presents the analysis of thermal pollution from a nuclear power plant located on the shores of the tidal South China Sea. The paper noted that it is more difficult to evaluate the effect of thermal pollution from a power plant in tidal seas than in non-tidal seas.

At present, there are five operating NPPs on the shores of the Baltic Sea: two Swedish (Forsmark NPP, electric capacity of 3210 MW; Oskarshamn, 2308 MW), two Finnish (Loviisa NPP, 1020 MW; Olkiluoto NPP, 1760 MW), and one Russian (Leningrad NPP, 4000 MW). Two of them (Forsmark and Olkiluoto) are located on the coast of the Bothnian Sea. On 19 January 2016, the construction of NPP “Hanhikivi-1” with the capacity of 1200 MW was started. This event had been preceded by examination of hydro-meteorological conditions in the area of construction, which included not only the estimation of extreme conditions in the vicinity of the Hanhikivi peninsula (Pyhäjoki municipality) in the Bothnian Bay in the Baltic Sea but also the possible impact of the future power plant on the marine environment in this area.

In particular, the evaluation of extreme weather and sea events in the Bothnian Bay in the case of the absence of the NPP was carried out by the Swedish Meteorological and Hydrological Institute (SMHI, 2014). To find the extreme values of sea temperature and water level, SMHI used long time series of observations of these characteristics, complemented by the results obtained by High Resolution Operational Model for the Baltic Sea (HIROMB) runs with a horizontal resolution of 1 $\times$ 1 nautical mile and a vertical resolution of 4–8 m. The analysis was limited by the period of 1990–2011 in the case of temperature and by the whole period of observations in the case of water level. A similar assessment was made by the company Fennovoima Ltd., making use of only available observational data in the area of the planned construction (Fennovoima report, 2015; Alenius, 2015). Both estimates are based on statistical methods by which the maximum and minimum values and the extreme probability and frequency were calculated. As for the EIA of the proposed NPP, those estimates were usually based on observational data, as well as on simple dispersion models for jet propagation in liquid and expert assessments (Fennovoima report, 2014); thus, both estimates can be considered as preliminary and indicative. The method proposed below, which is based on three-dimensional modeling of hydrodynamics on ultra-high-resolution grids, provides more reliable assessments.

The purpose of this study was twofold: (1) to estimate the possible extreme marine phenomena in this region (wind waves, sea level changes) and (2) to estimate the impact of the NPP on the local thermal regime in the future. To do this, we used different hydrodynamic models. Results of this examination are presented below.

Methods

When evaluating extreme hydrometeorological conditions in the study area of the sea, an empirical approach based on statistical analysis of the available time series of individual characteristics (wind speed, sea level, water temperature, ice cover characteristics, and others) is traditionally used. This approach has at least two obvious limitations: (1) the accuracy of such estimates is highly dependent on the duration of the observation period and (2) it cannot be used for areas where such data are not available. In this paper we propose a method for estimating extreme hydrological conditions in the study area based on mathematical modeling of the general ocean circulation. In the present study, the extreme hydrological conditions in the area of the future NPP “Hanhikivi-1” are estimated in the Bothnian Bay in the Baltic Sea. The general scheme of the calculation is as follows: (1) we perform a model run to simulate the general circulation for the entire Bothnian Bay on the coarse grid for the selected period (from 2010 to 2015, with double repetition of the year 2010 under one and the same external forcing) including the cold (2010) and warm (2014) years in the situation of the absence of the NPP. The model performance is verified through the quality of simulation of temperature and sea ice area fields. (2) We assess the extreme possible sea level in the NPP area on the basis of runs on the coarse grid. For this purpose the real situations of extreme storm surges during the selected time period were chosen (3) to assess the impact of the NPP on the temperature and sea ice fields around it on the basis of performing the calculations in case of the absence of the NPP (“background” scenario) and its presence (“predictive” scenario); runs are performed on the nested (fine) grid covering the neighborhood of the NPP for both cold and warm years. Real situations of extreme storm surges during the selected time period were also chosen (4) to assess the largest wind waves by performing the numerical experiments first on the coarse grid and then on the fine grid, with the prescribing boundary conditions from the solutions on the coarse grid.

The general scheme of numerical experiments is shown in Fig. 1. The models used and the details of the numerical experiments above are specified below.

The general scheme of the model calculations.

[Figure omitted. See PDF]

Bathymetry of the Bothnian Bay (a) and the area off the Hanhikivi peninsula (b). The location of the NPP, oceanographic stations (BO3, F3, PP2, PP4, PP5), coastal sites of sea level measurements (Ratan, Furuögrund, Kalix, Kemi, Raahe, Pietarsaari), and also the position of cross section (red line) and SWAN output point (W) are presented. The inset in fragment (a) shows the position of the Bothnian Bay in the Baltic Sea.

[Figure omitted. See PDF]

Circulation model

A three-dimensional numerical model based on the Princeton Ocean Model (POM) (Blumberg and Mellor, 1987; Mellor, 2004) was used to simulate the circulation pattern and thermal regime in the Bothnian Bay. It is a model with a $\mathit{\sigma }$ coordinate in the vertical direction that allows a smooth representation of the bathymetry and provides a better simulation of the currents in the bottom boundary layer as compared with some $z$ and isopycnal models. Also, such $\mathit{\sigma }$ models reproduce different thermodynamic effects caused by the nonlinearity of the equation of state rather well. That is why such models have been widely used for the simulations of coastal and estuarine dynamics. As for the turbulence closure parameterization, POM makes use of the well-tested and reliable schemes of Mellor–Yamada 2.5 (Mellor and Yamada, 1982) for the vertical and Smagorinsky (Smagorinsky et al., 1965) for the horizontal turbulent mixing. The abovementioned advantages of the model, along with the significant experience of the authors with its use in a number of successful applications for other basins in similar studies (Ryabchenko et al., 2008), were the main reasons why POM was chosen for the present study. The model allows the simulation of three-dimensional fields of current, water temperature and salinity, two-dimensional elevation field in an area of interest when external influence is specified. This influence includes wind forcing, sea surface atmospheric pressure, air temperature and humidity in the near-water layer, cloudiness, precipitation rate, and open boundary conditions.

Sea ice and snow model

To model ice and snow distribution in any area, an advanced sea ice model with several different categories of ice was used (Haapala et al., 2005; Ryabchenko et al., 2010). The model distinguishes the sea ice as two main types: deformed and non-deformed. The non-deformed ice is divided into several subcategories, while the deformed ice consists of only ridged and rafted ice. The rafted ice exists when the ice thickness is equal to or less than 17 cm. Otherwise, the ice is considered ridged. Evolution of each type is described by ice concentration and mass equations. The ice thickness of each category varies due to advection, deformation, and thermodynamic processes. It is assumed that fast ice exists in areas with depth less than some specified value. However, in the present study fast ice was not considered, and the model operated with seven different categories of sea ice.

The evolution of snow cover thickness can be described as the interaction of the following main mechanisms (Lepparanta, 1983): precipitation, surface melting, compaction, and the formation of slush, which further transforms into snow ice. In the present model only precipitation and surface melting were taken into account (Ryabchenko et al., 2010).

Wind wave model

The wind wave model SWAN (Simulating Waves Nearshore) (Booij et al., 1999; Ris et al., 1999) was used in the present study. SWAN is a third-generation spectral wave model specifically developed for wind wave simulation in shelf and shallow coastal areas with complex shoreline configuration. The model can take into account wind forcing, depth-induced wave breaking, refraction, diffraction, ambient currents and sea level oscillations, bottom friction, white capping, wave quadruplets and triads, wave-induced setup, presence of sub-grid obstacles, vegetation and bottom mud layer, and the turbulent viscosity. The core element of SWAN is the numerical and efficient solving of the spectral wave action balance equation, which includes source terms representing the effects of generation, dissipation, and nonlinear wave–wave interactions.

Statistical characteristics of air temperature, atmospheric pressure, and wind speed in the surface layer of the atmosphere between 2010 and 2014, calculated on the average daily data.

Years	Air temperature, ${}^{\circ }$C		Atmospheric pressure, mb		Wind velocity, m s${}^{-\mathrm{1}}$
	Observations	HIRLAM	Observations	HIRLAM	Observations	HIRLAM
Mean annual value
2010	2.1	2.0	1012.2	1012.6	5.5	4.7
2011	4.6	3.7	1009.0	1009.5	6.3	5.5
2012	3.1	3.2	1010.3	1010.7	5.8	6.1
2013	4.7	4.7	1010.2	1010.6	6.2	6.1
2014	5.2	5.1	1012.8	1013.0	6.1	6.1
Minimum value
2010	$-$27.4	$-$26.5	981.8	981.6	1.5	1.1
2011	$-$25.6	$-$26.9	969.3	970.4	1.2	0.9
2012	$-$24.7	$-$26.1	973.2	973.4	1.7	1.9
2013	$-$19.2	$-$20.1	974.9	976.1	2.0	1.0
2014	$-$20.0	$-$21.9	976.5	977.3	1.7	1.6
Maximum value
2010	26.0	24.7	1046.8	1046.3	11.6	13.1
2011	25.6	20.4	1037.4	1038.0	18.9	17.2
2012	20.1	20.1	1056.0	1055.7	13.2	14.6
2013	21.4	21.3	1033.6	1033.9	17.7	18.1
2014	25.8	24.7	1040.6	1040.1	14.7	15.6
Standard deviation
2010	11.4	12.0	11.3	11.2	2.2	2.1
2011	9.9	9.0	12.8	12.7	2.8	2.9
2012	9.7	9.9	13.0	12.9	2.3	2.5
2013	9.4	9.2	10.7	10.7	2.5	2.8
2014	8.8	9.1	11.8	11.6	2.7	2.8

SWAN was used in a nonstationary mode with a time step equal to 10 min, with maximal iterations at each time step set equal to 10. According to the recommendations presented in the official SWAN manual, the directional resolution was set equal to 10${}^{\circ }$, and minimal and maximal frequencies were set equal to 0.04 and 1.0 Hz.

Areas of interest and model domains

The present study focuses on the Bothnian Bay (a northern part of the Baltic Sea) and the relatively small area off the Hanhikivi peninsula located at the eastern coast of the Bothnian Bay (Fig. 2). The bathymetry data were collected from different marine navigational charts, from the Baltic Sea Bathymetry Database (Baltic Sea Hydrographic Commission, 2013), and were provided from field observations. These data were linearly interpolated into grid nodes, using a low-frequency filter.

The Bothnian Bay model grid consists of 107 $\times$ 175 nodes horizontally and 25 $\mathit{\sigma }$ levels vertically, horizontal resolution being one nautical mile and minimum depth being 3 m. The model time step was set equal to 300 s.

The Hanhikivi area model grid is much more curvilinear in its shape (see Fig. 2b), consisting of 142 $\times$ 193 nodes horizontally and 12 $\mathit{\sigma }$ levels vertically. Its horizontal resolution is rather variable, with the minimum and maximum horizontal resolution being 35 and 180 m, respectively. The minimum depth was set equal to 1.5 m and the model time step was 0.5 s. The Hanhikivi model grid covers an area with a radius of approximately 9000 m from the future NPP Hanhikivi-1.

In addition to these abovementioned main model domains and their grids, we also built and used another grid, which covered both the Bothnian Sea and the Bothnian Bay, in order to assess the impact of incoming wind waves traveling from the Bothnian Sea northward into the Bothnian Bay through the open south boundary. This model grid (not shown) consisted of 101 $\times$ 250 nodes. The minimum depth was limited by the 3 m isobath, with a horizontal resolution of 2.2 km minimum and 5.3 km maximum. The main outcome obtained from the numerical experiments performed on this grid was that the impact of the incoming wind waves plays an important role only in the southern part of the Bothnian Bay and that its contribution to the wind wave pattern in the vicinity of the Hanhikivi peninsula is negligible. Such a result allowed us to concentrate only on the modeling efforts inside the Bothnian Bay without any need to include the Bothnian Sea in the model domain in all subsequent model runs, thus enhancing the model resolution and reducing the model total integration time.

Atmospheric forcing, boundary and initial conditions

Atmospheric forcing included air temperature and humidity, wind speed and direction, and cloudiness, which were obtained from the results of the atmospheric model HIRLAM (High Resolution Limited Area Model) (http://hirlam.org), provided to Russian State Hydrometeorological University by the Danish Meteorological Institute. HIRLAM's horizontal resolution is 11 km and the model's domain covers the North Atlantic and northern European regions. To compare HIRLAM results with measurements, observational data from the meteorological station Raahe were used. These data represent the 3 h monitoring data on air temperature and atmospheric pressure for the period from 1 January 2010 to 31 December 2014 and on the characteristics of wind speed for the period from 22 October 2010 to 31 December 2014. For the same period similar HIRLAM characteristics were made. The model and the observed values were averaged with a daily period. Table 1 shows the statistical characteristics of the observed and model values from year to year, and Fig. 3 shows a comparison of calculated and observed mean monthly values of meteorological characteristics in 2010–2014. As can be seen, there is a good agreement between the model and the observed values and thus the HIRLAM atmospheric forcing can be used in the simulation of the hydrodynamic regime of the Bothnian Bay near the future NPP location.

Temporal variability of mean monthly values of (a) air temperature (${}^{\circ }$C), (b) atmospheric pressure (mb), and (c) wind speed (m s${}^{-\mathrm{1}}$) in the near-surface layer of the atmosphere. HIRLAM (dashed curves) and observations (green curves).

[Figure omitted. See PDF]

The circulation model used in the present study calculates the momentum, heat, and salt fluxes at the air–sea and air–ice interfaces. The momentum fluxes are calculated traditionally as a quadratic friction law, making use of different drag coefficients for air–water, air–ice, and ice–water boundaries. The heat and salt fluxes are parameterized by taking into account the diurnal cycle of shortwave solar radiation (Parkinson and Washington, 1979; Ryabchenko et al., 2010).

At the open boundaries the coupled circulation and sea ice model assimilates the sea level, current velocity, water temperature and salinity, sea ice thickness and concentration, and snow thickness obtained from the results of HIROMB (Funkquist, 2001; http://www.smhi.se/en/research/research-departments/oceanography), a high-quality prognostic 3-D ocean circulation model recommended for all states in the Baltic Sea region by the Helsinki Commission (HELCOM). HIROMB covers the North Sea and the Baltic Sea, with horizontal resolution in the latter being 1 nautical mile and vertical resolution being 4 m (24 vertical levels). HIROMB data also include river runoff in the Baltic Sea region.

At the solid lateral boundaries (coasts) a no-slip condition and zero fluxes are specified for horizontal velocity, heat, and salt.

At the bottom the vertical component of current velocity, heat, and salt fluxes is set to zero. Bed shear stress is parameterized as a function of horizontal velocity at the near-bed model $\mathit{\sigma }$ level.

All the hydrological boundary conditions described above were applied for the Bothnian Bay circulation model, while for the Hanhikivi area circulation model the results obtained from the previous Bothnian Bay model runs were used as the boundary conditions, implementing a nesting technique.

It should also be noted that the hydrological regime of the area located relatively close to the Hanhikivi peninsula is considerably affected by the river Pohjoishaara situated nearby, with an annual runoff equal to 103 $\times$ 10${}^{\mathrm{7}}$ m${}^{\mathrm{3}}$. In the present study, for the Hanhikivi model the river Pohjoishaara runoff was distributed into 12 months in accordance with the known monthly averages of the nearest large rivers' runoff. The Pohjoishaara's water temperature was set to the nearest coastal area water temperature, while salinity was set to zero.

The initial conditions for the coupled circulation and sea ice model for the entire Bothnian Bay domain included water temperature, salinity, and sea level fields obtained from the HIROMB model.

Generally speaking, extreme sea levels in the Bothnian Bay are caused by storm winds, long waves, tides, low atmospheric pressure, seiches, and sea level rise of the world ocean. Tidal level oscillations in the Gulf of Bothnia are negligible (Lepparanta and Myrberg, 2009) and can be omitted in model simulations. The influence of moving centers of low atmospheric pressure has not been investigated in the present study. Still it can be assumed that their impact commonly appears jointly with the wind impact (SMHI, 2014). In order to simulate extreme sea level oscillations in the vicinity of the NPP Hanhikivi-1, we considered the situations with constant (both in speed and direction) maximal possible wind blowing long enough to establish an equilibrium state and under the influence of sea level change caused by the long waves coming from the Baltic Sea. As indicated above, such a simplistic approach to the evaluation of extreme sea level in the area of interest gives results in good agreement with the estimates of extreme values of sea level according to the observations over almost 100 years.

For the SWAN model the initial condition of no waves in the entire domain was adopted. At the solid boundaries the model sets the condition of full wave energy absorption and no wave energy generation. At the open boundaries in the case of outgoing waves, they can leave the area freely (radiation condition), while in the case of incoming waves there can be two options: (1) either no waves come into the model domain (an option used for the Bothnian Bay model) or (2) incoming wave spectrum is specified along the open boundary, which has been obtained from the previous Bothnian Bay simulations (an option used for the Hanhikivi model).

The model estimates of wind wave characteristics took into account the possible presence of sea ice cover in the Bothnian Bay, which in reality hinders or prevents the free generation and propagation of surface waves and limits the wave fetch if some part of the basin is covered with ice. The model run for 2014 was fulfilled for the entire year, which is why the inclusion of sea ice in the wind wave model was necessary. We assumed the isoline of ice concentration equal to 0.5 as an edge of the ice cover, which was in fact some kind of simplification but still appeared to be an effective way to limit the wave fetch in the presence of sea ice in the Bothnian Bay. All required data for the SWAN model (sea level oscillations, currents, sea ice) had been calculated in advance by the coupled circulation model before being used in the SWAN model.

Comparison of POM- and HIROMB-calculated sea levels at Raahe station located near the NPP Hanhikivi-1 for the storm surge periods of 14–16 October 2010 (a) and 5–7 December 2015 (b).

[Figure omitted. See PDF]

Comparison of the computed temperature profiles by POM and HIROMB with observations provided by BED in the Bothnian Bay. Locations of the stations are presented in Fig. 2a. (a) Station BO3, 18 January 2010; (b) station F3, 19 January 2010; (c) station BO3, 6 June 2010; (d) station F3, 30 June 2010.

[Figure omitted. See PDF]

Verification of the models

Taking into account the main objective of this study – to estimate extreme values of the sea level and height of wind waves in the vicinity of the NPP Hanhikivi-1 and estimate the maximum thermal pollution produced by the NPP – the proposed models were verified with respect to the sea level, significant wave height, sea water temperature, and sea ice area against all available observational data for the selected period 2010–2015. These data include data on sea level at six stations (Ratan, Furuögrund, Kalix, Kemi, Raahe, Pietarsaari) on the shores of the Bothnian Bay (http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=INSITU_BAL_NRT_OBSERVATIONS_013_032); vertical temperature profiles for hydrographic stations BO3, F3 (data from BED http://www.balticnest.org/bed), and PP5 (Fennovoima report, 2014); data on the significant wave height at stations PP2 and PP4 (Fennovoima report, 2013); and the average monthly data on the area of the ice cover (http://www.aari.ru/).

Circulation model verification

Sea level. Model verification with respect to the sea level was carried out using the series of hourly sea level values for the period from 00 h of 1 January 2010 to 23 h of 31 December 2015, obtained from observational data and modeling results at the six stations mentioned above. Statistical characteristics (correlation coefficient $R$ and standard deviation $\mathit{\sigma })$ of sea level at the coastal stations (Table 2) indicate that both models reproduce the level with high enough accuracy everywhere in the Bothnian Bay. The observed and calculated sea levels for the short periods of 14–16 October 2010 and 5–7 December 2015 at the station Raahe nearest to the NPP are presented as an example (Fig. 4). Although the two models underestimate the amplitude of sea level changes, in general they give consistent results and reproduce the sea level changes with sufficient accuracy. These periods were chosen as periods of strong storm surges causing large fall and rise of the sea level at the Raahe station. The extreme high and low sea level values are estimated in Sect. 3.1.2 using information about the initial conditions for the hypothetical calculations of extreme values of levels just for the beginning of these two periods.

Statistical characteristics (correlation coefficient $R$ and standard deviation $\mathit{\sigma }$) of water level at coastal stations in the Bothnian Bay in 2010–2015 according to observations and results from HIROMB and POM (this study).

Characteristic	$R$,	$R$,	$\mathit{\sigma }$, m	$\mathit{\sigma }$, m	$\mathit{\sigma }$, m
	data – POM	data – HIROMB	data	POM	HIROMB
Station
Ratan	0.88	0.88	0.24	0.21	0.22
Furuögrund	0.88	0.89	0.25	0.22	0.24
Kalix	0.88	0.88	0.27	0.24	0.25
Kemi	0.89	0.90	0.28	0.25	0.27
Raahe	0.88	0.89	0.23	0.26	0.25
Pietarsaari	0.80	0.82	0.21	0.22	0.23

Temperature. Comparison of the results of temperature simulation on the coarse grid with the available data of field observations in the Bothnian Bay from BED were performed for two stations: BO3 (64.30${}^{\circ }$ N, 22.35${}^{\circ }$ E) and F3 (65.17${}^{\circ }$ N, 23.23${}^{\circ }$ E). At these stations for the considered period in open access there are 40 vertical profiles of temperature: 12 for BO3 and 28 for F3. Comparison of mean values and standard deviations for the observed vertical profiles of temperature with similar characteristics obtained on a 2 nm grid using HIROMB and on a 1 nm grid using POM (this work) is presented in Table 3. Figure 5 shows the examples of profiles for both stations in January and June 2010.

Statistical characteristics (mean value m and standard deviation $\mathit{\sigma }$) of vertical temperature profiles at the oceanographic stations BO3 and F3 in the Bothnian Bay in 2010–2012 according to BED observations, results from HIROMB, and POM (this study).

Station	Station	Station	$m$, ${}^{\circ }$C	$m$, ${}^{\circ }$C	$m$, ${}^{\circ }$C	$\mathit{\sigma }$, ${}^{\circ }$C	$\mathit{\sigma }$, ${}^{\circ }$C	$\mathit{\sigma }$, ${}^{\circ }$C
number	name	date	data	HIROMB	POM	data	HIROMB	POM
1	BO3	2010-01-18	1.29	2.90	2.02	1.36	2.64	1.77
2	BO3	2010-01-18	0.60	1.44	0.53	0.56	1.43	0.57
3	BO3	2010-06-06	1.73	2.20	1.89	1.10	1.14	1.11
4	BO3	2010-08-22	4.05	6.50	7.11	4.80	5.16	4.56
5	BO3	2011-01-24	0.82	1.47	0.76	0.88	1.10	0.67
6	BO3	2011-05-29	1.30	1.92	1.76	1.01	1.56	1.23
7	BO3	2011-08-21	5.76	6.37	8.34	4.65	6.01	4.64
8	BO3	2012-05-31	1.99	2.71	2.93	1.07	1.49	1.17
9	BO3	2012-08-23	6.13	7.83	8.38	4.13	4.85	3.98
10	BO3	2010-12-02	3.61	2.91	3.01	0.44	1.14	1.32
11	BO3	2011-12-11	3.30	2.92	3.53	0.52	0.80	0.17
12	BO3	2012-12-10	3.58	3.93	4.08	0.39	1.21	1.21
13	F3	2010-01-19	1.74	2.45	2.41	1.53	2.22	2.47
14	F3	2011-01-24	0.95	1.14	1.27	0.93	1.08	1.29
15	F3	2010-03-05	0.41	2.14	1.81	1.12	2.28	1.79
16	F3	2010-03-23	0.73	1.91	1.92	1.09	2.04	1.95
17	F3	2010-05-05	0.92	2.44	2.29	1.11	1.94	1.67
18	F3	2010-05-19	1.29	2.22	2.20	0.88	1.78	1.44
19	F3	2010-06-09	3.74	3.53	2.94	2.79	0.34	0.72
20	F3	2010-06-30	5.30	5.76	5.07	4.24	2.71	2.75
21	F3	2010-08-04	7.05	7.82	7.63	6.43	5.33	4.66
22	F3	2010-09-01	6.49	6.43	7.67	6.02	5.16	4.29
23	F3	2010-10-27	4.01	5.39	6.18	1.30	1.26	0.34
24	F3	2010-12-08	2.82	2.28	2.84	1.09	0.61	1.05
25	F3	2011-03-15	0.15	1.44	1.24	0.50	1.13	1.02
26	F3	2011-05-18	1.69	1.03	1.13	1.15	1.13	0.76
27	F3	2011-06-08	3.66	3.11	2.60	2.03	1.57	1.17
28	F3	2011-07-06	6.41	5.58	5.36	5.50	4.03	3.66
29	F3	2011-07-20	7.25	6.98	7.13	6.59	5.41	5.00
30	F3	2011-08-03	9.02	8.00	8.37	8.49	6.25	5.69
31	F3	2011-08-17	6.41	7.03	7.35	5.50	5.48	5.03
32	F3	2011-09-07	9.02	6.52	7.44	8.49	4.94	4.52
33	F3	2011-10-26	5.42	5.53	6.89	2.70	2.18	0.84
34	F3	2011-12-07	4.62	4.56	4.89	0.44	0.22	0.23
35	F3	2012-01-18	2.22	1.88	2.46	0.12	0.12	0.52
36	F3	2012-03-20	2.22	$-$0.05	0.16	0.12	0.04	0.24
37	F3	2012-04-25	0.27	0.29	0.20	0.25	0.38	0.20
38	F3	2012-05-09	0.67	0.53	0.55	0.20	0.21	0.10
39	F3	2012-08-29	7.67	6.60	7.10	6.66	5.68	4.99
40	F3	2012-12-22	1.73	1.99	2.51	2.02	1.35	1.55

Analysis of Table 3 shows that when considering all profiles, with respect to average temperature values $m$, POM is closer to the observation data than HIROMB in 58 % of cases, and with respect to the standard deviation $\mathit{\sigma }$, POM is closer only 50 % of cases. When considering the stations separately, for station BO3 POM is better with respect to $m$ and $\mathit{\sigma }$ in 58 and 75 % of cases, while the estimates for station F3 are 57 and 39 %, respectively. In general, from this comparison it can be concluded that the performance in reproducing the water temperature in the Bothnian Bay on coarse grids for both models is about the same.

Comparison of the computed temperature profiles by POM on coarse and fine grids (blue and red curves, respectively) with observations (black curves) (Fennovoima report, 2014) at station PP5 near the NPP Hanhikivi-1 in the Bothnian Bay for (a) 8 June 2010; (b) 3 August 2010. Location of station PP5 is presented in Fig. 2b.

[Figure omitted. See PDF]

Obviously, the coarse grid with a step of 1–2 nm is not able to reproduce the spread of plumes of warm water from the NPP Hanhikivi-1 discharge point because the characteristic length scale of the plume is of the order of 1 km. Proposed in this study, the fine nested grid for the vicinity of NPP Hanhikivi-1 with steps from 35 to 180 m will solve this problem, i.e., will accurately reproduce the size and shape of thermal pollution plumes around the station. It is important that the solution on the fine grid is significantly closer to observational data in comparison with the solution on the coarse grid (see Fig. 6).

Sea ice model verification

The comparison of computed sea ice thickness and compactness showed that the model results were in general accordance with observed values (http://www.aari.ru/), though some overestimation in ice thickness can be observed (Fig. 7).

Modeled (a) and observed (b) distribution of sea ice thickness in the Bothnian Bay during the period 27 February–1 March 2011. Ice map in panel (b) is from the Arctic and Antarctic Research Institute Center “Sever” (St. Petersburg).

[Figure omitted. See PDF]

Wind wave model verification

Verification of the results of the SWAN model (in terms of significant wave height, SWH) against observational data (Fennovoima report, 2013) for two points (PP2 and PP4) located in the vicinity of the Hanhikivi peninsula showed that in general the model correctly simulated wind wave characteristics (Fig. 8). The main discrepancies could be caused by inaccuracy in dealing with ice cover in the wave model and/or ice cover modeling itself.

Comparison of SWH observations and model results for the periods (a) 6 November 2012–17 October 2013 and (b) 30 August 2012–17 October 2013.

[Figure omitted. See PDF]

(a) Bothnian Bay imaginary axis (red line), wind direction in numerical experiments (blue arrows), and Hanhikivi peninsula location (red square); (b) dependence of modeled SWH near the NPP upon wind direction; (c) dependence of modeled SWH near the NPP upon wind speed for the most dangerous wind direction (310${}^{\circ }$).

[Figure omitted. See PDF]

Summarizing, we can say that the proposed POM-based modeling system reproduces the principal characteristics of hydrodynamic regime (level, water temperature, altitude wind waves, sea ice) on the coarse grid at least no worse than the best model of the Baltic Sea HIROMB and gives a considerably better description of the temperature field on the fine grid. An advantage of POM that is important for the prognostic runs is the fact that the POM, unlike HIROMB, does not assimilate observational data.

Results

Modeling of storm events in the vicinity of the future NPP

Wind waves

With an aim to assess the maximal possible SWHs in the area off the Hanhikivi peninsula and to investigate the main features of wind wave fields in different hydrometeorological situations, a number of numerical experiments with different analytically specified wind speeds and directions was carried out. The area of interest is located along the eastern coast in the Bothnian Bay; thus, we expected that the maximal wind waves would be generated by western winds. If one drew a hypothetical line along the Bothnian Bay, then it would have the direction from southwest to northeast with the angle equal to approximately 50${}^{\circ }$ relative to the geographical west–east parallel. Using such a line as a coordinate axis, we performed seven main numerical experiments (A–G), making use of the SWAN model for wind directions varying from 0 to 180${}^{\circ }$ relative to this basic line, with angle discretization equal to 30${}^{\circ }$ and wind speed constant and equal to 10 m s${}^{-\mathrm{1}}$ (Fig. 9a).

These wind directions cover the entire range of possible directions capable of producing high waves in the vicinity of the future NPP Hanhikivi-1. All other wind directions outside the abovementioned range will lead to the occurrence of wave shadow zone near the Hanhikivi peninsula. Wind duration was 24 h. After the most dangerous wind direction had been determined, we used that direction and varied wind speeds, holding the direction constant. For that experiment the wind duration was also 24 h.

Model calculations of wind waves have shown that the most dangerous in terms of the generation of wind waves in the NPP area are the west and northwest wind, with the directions of 280 and 310${}^{\circ }$ (experiments E and D, respectively; Fig. 9b). Maximal SWH in the Bothnian Bay near the NPP for that wind direction at a wind speed of 10 m s${}^{-\mathrm{1}}$ is about 1.2–1.4 m after 24 h of wind. Changing the wind speed for the determined most dangerous wind direction (310${}^{\circ }$) allowed the assessment of the values of the highest possible wind waves near the NPP (Fig. 9c).

The model results allowed the estimation of the values of SWH in both the entire Bothnian Bay and in the small area near Hanhikivi peninsula during different external wind forcing. SWH for the most dangerous wind direction and wind speed of 15 m s${}^{-\mathrm{1}}$ was 2.5–3.0 m and increased up to 7.0 m in the open part of the Bothnian Bay for the wind speed of 30 m s${}^{-\mathrm{1}}$ (not shown). Nevertheless, in the very vicinity of the coast near the Hanhikivi peninsula, the SWH decreases dramatically due to relatively shallow depths (see Fig. 2b). It is also interesting to note that a considerable increase in SWH in this area in the model experiments appeared during the increase in wind speed from 10 to 15 m s${}^{-\mathrm{1}}$, while further wind speed increase did not show such rapid SWH growth (Fig. 9c).

In addition to the model experiments with theoretical atmospheric forcing, a model run for the real meteorological forcing of 2014 was carried out for the whole year. Just as in all other numerical experiments, both models (Bothnian Bay model and Hanhikivi area model) were used: the Bothnian Bay model produced all required wind wave characteristics to be used as boundary conditions in the Hanhikivi area model. The results obtained for 2014 show that the highest waves in the Bothnian Bay occurred during the autumn period, with SWH reaching 4.0 m. During the winter period of 2014, SWH reached 1.5–2.0 m at ice-free areas of the Bothnian Bay. However, most of the time during moderate wind speeds, SWH was at most 0.5–1.5 m.

A situation with rather large simulated SWH occurred on 27 September 2014 and is presented in Fig. 10 in order to show in detail the spatial distribution of wave heights in this area and the influence of bathymetry upon them. The region of the reduction of SWH due to wave breaking and bottom friction is clearly visible and in general coincides with the 5 m isobath. It should be noted that we used the parameterization for bottom friction implemented in the SWAN model, which takes into account the size of bed ripples calculated by making use of available field measurements of grain sizes.

Modeled SWH distribution on 27 September 2014.

[Figure omitted. See PDF]

Wind waves in the area near the Hanhikivi peninsula can be characterized by the time series of computed SWH at the location marked by the letter “W” in Fig. 2b. The time series of SWH in the Hanhikivi area for the whole of 2014 is presented in Fig. 11.

Modeled SWH time series in the Hanhikivi area at the location mark with “W” in Fig. 2b.

[Figure omitted. See PDF]

(a) Sea level rise in the Bothnian Bay during the model run with the wind speed 30.2 m s${}^{-\mathrm{1}}$ and a high level of the storm surge; (b) temporal evolution of sea level rise near the NPP area. $x$ axis corresponds to the time in hours from the start of the wind.

[Figure omitted. See PDF]

(a) Sea level rise in the Bothnian Bay during the model run with the wind speed 30.2 m s${}^{-\mathrm{1}}$ and low level of the storm surge; (b) temporal evolution of sea level rise near the NPP area. $x$ axis corresponds to the time in hours from the start of the wind.

[Figure omitted. See PDF]

Extreme sea levels

To calculate the extreme sea levels in the vicinity of the NPP Hanhikivi-1, the maximum possible wind velocity was set equal to 30.2 m s${}^{-\mathrm{1}}$ (having 0.01 exceedance probability), defined according to the observations in the NPP area (Fennovoima Report, 2015). The wind velocity field calculated using the HIRLAM model was enhanced to reach this maximum wind velocity without any changes in wind direction. The periods of 5–7 December 2015 and 14–16 October 2010 were chosen to simulate the storm surges, causing the extreme high and low sea level rise, respectively, because during those periods the real extreme storm surges in the Bothnian Bay were observed. The years 2010–2015 were chosen because we had data for the verification of the model mainly from that period. The results for the entire Bothnian Bay are presented in Fig. 12 for the high level of the storm surge and in Fig. 13 for the low level of the storm surge. All computations were performed by making use of the Bothnian Bay coupled circulation model only, without using the nesting technique and the Hanhikivi area model.

Ice thickness distribution (monthly mean) in the vicinity of the NPP Hanhikivi-1 in February for conditions of the warm year 2014. (a) Background scenario; (b) predictive scenario.

[Figure omitted. See PDF]

Ice thickness distribution (monthly mean) in the vicinity of the NPP Hanhikivi-1 in February for conditions of the cold year 2010. (a) Background scenario; (b) predictive scenario.

[Figure omitted. See PDF]

According to the simulations of possible sea level changes in the Bothnian Bay, extreme sea level values in the vicinity of the future NPP at a constant wind of 30.2 m s${}^{-\mathrm{1}}$ had a maximum of 248 cm and minimum of $-$151 cm (see Figs. 12b and 13b). This is in a good agreement with the sea level data for the period 1922–2015 at the Raahe station (Finland) according to the estimates of SMHI (SMHI, 2014) in which the maximal and minimal sea level there could reach up to 250.4 and $-$150 cm, respectively, once in 10${}^{\mathrm{8}}$ years.

Assessment of thermal pollution of the water by the NPP

To assess the possible impacts of the NPP Hanhikivi-1 on the local thermal regime, two scenario runs were performed:

• A background scenario simulating the natural conditions in the absence of the NPP was performed.

• A predictive scenario was performed. The NPP has been built and is operating with the temperature of heated discharged water set equal to 12 ${}^{\circ }$C above the ambient water temperature at the point of water discharge, the discharge being 45 m${}^{\mathrm{3}}$s${}^{-\mathrm{1}}$.

Runs were performed for a cold year (2010) and a warm year (2014). The atmospheric characteristics necessary for calculating the fluxes of moment, heat, and moisture at the air–water boundary were set according to the atmospheric HIRLAM model, with a time resolution of 1 h. To set the boundary conditions on the open boundary of the Bothnian Bay, we used the data from the HIROMB model (sea level, water and ice current velocities, temperature, salinity, ice thickness and compactness, snow thickness). The average discharge of the Pohjoishaara River was set equal to 33 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$.

The vertical structure of the temperature field along the section from the discharge point to the north for 4 April 2010: (a) background scenario; (b) predictive scenario. Temperature scales for (a) and (b) are very different.

[Figure omitted. See PDF]

In natural conditions, the water of the gulf around the Hanhikivi peninsula was covered with ice from the beginning of December to the beginning of May in 2010 (the cold year) and from the beginning of January to the beginning of April in 2014 (the warm year). The highest background temperature in the cold and warm years was achieved in July and August, respectively. The thermal regime of the basin in the vicinity of the points of water intake and water discharge is almost identical. In general, the spatial variations of mean monthly temperature in the area of the Hanhikivi peninsula limited by a radius of 2 km are small, not exceeding 0.6 ${}^{\circ }$C for the sea surface temperature (SST) and 1.2 ${}^{\circ }$C for the temperature of the deep layer. The main difference between the cold 2010 year and the warm 2014 year is a longer winter period in 2010.

Permanent discharge of warm water in the case of the operating NPP will lead to a permanent polynya near the northern tip of the Hanhikivi peninsula, resembling an ellipse with axes 1.5 $\times$ 6 km in warm winter conditions, which stretches to the north (Fig. 14b). In cold winter conditions the polynya has a rounded shape, the boundary of which is located at 2–3 km from the water discharge point (Fig. 15b).

The difference between background and predictive model runs is clear. A thermal plume (plume of heated water) emerges in the area of water discharge. Its spatial expansion and propagation mainly depends upon the wind speed and direction above the ice-free water surface and upon the current velocity and direction during ice cover periods. Figures 14–15 demonstrate the influence of the heated water discharge upon the ice cover distribution and thickness for both warm (2014) and cold (2010) years, respectively.

A vertical structure of water both in natural conditions and after the construction of the NPP was also investigated. A vertical structure changes significantly when large amounts of heated water have been discharged. Figure 16 shows the example of vertical cross sections of the temperature field calculated for the cold year conditions. The position of the cross section is presented in Fig. 2b. It starts from the heated-water discharge point and stretches to the north. The difference in depth for this one and the same profile is due to the changes in bathymetry caused by the planning hydrotechnical works near the future station.

The thermal regime in the vicinity of the water discharge point becomes completely different: SST deviations from background values are maximal in the 0–250 m zone, where they reach approximately 10 ${}^{\circ }$C in winter and spring and 8 ${}^{\circ }$C in summer and early autumn in the cold year (2010). SST deviations decrease with the increasing distance from the water discharge point and reach the minimum value of 0.3 ${}^{\circ }$C in the zone of 1500–2000 m. The bottom-water temperature is also maximal in the 0–250 m zone, where its deviations amount to 5 ${}^{\circ }$C in winter and spring and 3 ${}^{\circ }$C in summer and early autumn. In the cold year conditions the thermal regime in the case of the operating NPP will be changed qualitatively as well as in warm year, but these changes (compared to the background scenario) in the vicinity of both the intake point and the discharge point will be stronger.

Discussion

Assessment of the scale of the thermal effects that could arise due to the influence of the NPP Hanhikivi-1 on the local thermal regime has been obtained for anomalously warm and cold years. These years for the Hanhikivi area were identified as a result of the statistical analysis of long-term variability of meteorological parameters (with 3 h resolution) for the period from 1993 to 2014, observed at the meteorological station Raahe Lapaluoto. The data were provided by Fennovoima (Fennovoima report, 2015). It was found that the coldest winter for the period was observed in 2010. The hottest summer for the period from 1993 to 2014 was observed in 2002 and for the period from 2004 to 2014 the hottest summer was observed in 2014. Choosing 2014 as the abnormally warm year was dictated by the availability of observational data needed for verification of the models used: the number of available data for 2014 were greater than for 2002. If 2002 were considered as an abnormally warm year, the difference between the cold and warm years would be more. In any case, the resulting estimates should be viewed only as typical ones that determine the order of magnitude.

Another restriction of this study is connected with prescribing the constant temperature difference (12 ${}^{\circ }$C) between the discharged cooling water and the water temperature in the bay. This assumption, not taking into account seasonal variations in the temperature of water environment, can significantly affect the final result. For example, at the Beloyarsk NPP operating in the Ural region in the Russian Federation, this difference varies between 9 ${}^{\circ }$C in summer and 12–13 ${}^{\circ }$C in winter (Vereshchagina et al., 2013). Finally, this temperature difference can be set when the plant is operating.

The main missing physical mechanism in the coupling of the models used in the present study is the influence of the wave bottom boundary layer upon the current bottom boundary layer. The transfer of SWAN results into the POM model was not made during the computations, only from POM to SWAN. However, from our point of view, this interaction between bottom boundary layers has less impact on the final solution than the influence of meteorological forcing and model estimates of other hydrological parameters such as sea level, sea ice distribution, and currents.

Concluding remarks

• The approach used in this study based on the method of nested grids is well-known in oceanography. Here it has been used together with consideration of (1) model situations (setting the maximum possible wind of a certain direction in the selected periods of storm surges, prescribing direction and wind speed in the simulation of wind waves) aimed at evaluating the extreme changes in sea level and wind wave parameters in the selected area – the construction area of the NPP Hanhikivi-1 in the Bothnian Bay, Baltic Sea – and (2) the scenarios of warm and cold years for the detailed assessment of the thermal pollution of the NPP's neighborhood. One important feature of the nested grids used should be emphasized: the grid vertical structure does not change when going from the coarse to the fine grid. This avoids the situation with arising unstable stratification in fields interpolated on the fine grid and ensures the absence of numerical noise, which often causes the instability of computing. The scientific value of this approach is the fact that, unlike traditional statistical estimates of the extreme values of marine characteristics and their repeatability by observations from meteorological stations, it can be used in the local areas where the duration of the time series of observations is small (for example, in high-latitude Arctic seas).

• The most important results from the engineering point of view for the neighborhood of the NPP Hanhikivi-1 are as follows. Model calculations of wind waves have shown that the most dangerous (in terms of the generation of wind waves in the NPP area) is a northwest wind with the direction of 310${}^{\circ }$. The maximum height of the waves in the Bothnian Bay near the NPP for this wind direction with wind velocity of 10 m s${}^{-\mathrm{1}}$ is 1.2–1.4 m. According to the model estimates, the highest possible level of the sea near the NPP is 248 cm, and the minimum level is $-$151 cm, for the western and eastern winds, respectively. An important feature of thermal pollution around the station is that at the point of water intake for cooling the reactor, in the warm year the sea surface temperature in winter is not lower than the water freezing point and the ice cover is not formed, while in the cold year in the period from the end of December to March in this area there is ice cover with thickness of a few tens of centimeters. Thus, the extreme values of significant wave height, sea level, and water temperature in the vicinity of the future NPP Hanhikivi-1 can be significant and sometimes catastrophic.

• Numerical experiments for the cold (2010) and the warm years (2014) showed that the permanent release of heat into the marine environment from the operating NPP for the cold year will increase the temperature in the upper layer in the 250 m zone (from the heated water discharge point) by 10 ${}^{\circ }$C (in average) in winter and spring and by 8 ${}^{\circ }$C in summer and early autumn. In the bottom layer of the same zone, it will increase the temperature by 5 ${}^{\circ }$C in winter and spring and 3 ${}^{\circ }$C in summer and early autumn. For the warm year, these temperature changes will be smaller. In both cases, in the 1 km vicinity of the station throughout the whole year there is a vertical temperature structure with a pronounced thermocline at a depth of 1–2 m. Ice cover in both cases will disappear in a 2 km vicinity of the NPP. The warming of the marine environment is important to assess possible changes in the functioning of the natural aquatic ecosystems, including fish, when commissioning the NPP, as well as in choosing a site for the creation of aquaculture farms.

• According to the estimates, the scale of the thermal impact of the NPP on the local thermal regime is substantial and therefore this impact should be taken into account when assessing local climate changes and marine environmental impact in the future.

The data presented in this study are available in NetCDF format (https://swiftbrowser.dkrz.de/public/dkrz_035d8f6ff058403bb42f8302e6badfbc/Bothnia/, Dvornikov et al., 2017).

The authors declare that they have no conflict of interest.

Acknowledgements

This study was supported by grant 14-50-00095 of the Russian Science Foundation, by the grant 16-55-76021 of the Russian Foundation for Basic Research and by the German Federal Ministry of Education and Research (BMBF) under research grant 01DJ16016. We would like to thank the Finnish Meteorological Institute for kindly providing us with meteorological and sea level data at the Raahe station. We also thank the two anonymous reviewers for their helpful comments and suggestions. Edited by: A. Rutgersson Reviewed by: two anonymous referees