Assessment of Tidal Current Energy Resources in the Pearl River Estuary Using a Numerical Method

Abstract

In this paper, a numerical method is employed to assess tidal current energy resources in the Pearl River Estuary, China. The numerical model for tidal current simulation in the estuary is developed based on the MIKE 21 model, which enables numerical simulations in estuaries, coastal areas, and oceans. The model has a grid resolution that varies from about 2500 m at the open boundary to 500–1000 m inside the estuary. Extensive model validation is performed by comparing the model predictions with field observations of tidal level and velocity at various stations in the Pearl River Estuary. The tidal characteristics are thoroughly analyzed. Energy fluxes and power densities are calculated along selected cross sections to evaluate the feasibility of tidal energy development in the Pearl River Estuary. The results indicate that the distribution of annual average tidal current power density in the Pearl River Estuary generally aligns with the spatial distribution of tidal currents. The annual average power density of tidal energy is typically below 0.10 kW/m². The theoretical potential of tidal current energy resources in the Pearl River Estuary is assessed to be approximately 11,000 kW.

1. Introduction

The atmospheric concentration of carbon dioxide reached a new record of 4.132 × 10⁻⁴ mg/L in the 2021 Greenhouse Gas Bulletin published by the World Meteorological Organization (WMO), which corresponds to an increase of 149% compared to the pre-industrial level in 1750 [1]. Reducing global carbon emissions and addressing the climate and environmental issues caused by global warming has become major challenge faced by all of humanity. Due to its abundant reserves, ocean renewable energy has become an important means to address this challenge. Tidal current energy, which is generated by seawater’s periodic horizontal motion caused by the moon’s gravitational forces and the sun, is an important component of ocean renewable energy. Compared to other forms of marine energy, tidal energy has strong regularity and predictability (Lamy and Azevedo, 2018) [2]. Tidal energy conversion devices are typically installed on the seabed or float on the ocean surface, minimizing environmental impacts on ocean ecosystems and requiring minimal land resources. Tidal energy has a higher energy density compared to wind energy (approximately four times) and solar energy (approximately thirty times) [3]. The utilization of tidal current energy resources is a significant focus of energy research worldwide.

The characteristics of tidal currents and the tidal current energy resources are the foundation for designing and improving tidal current energy converters, as well as important considerations for tidal current energy development and site selection. To harness tidal current energy for electricity generation, it is essential to understand the tidal characteristics and the level of tidal current energy resources in the potential installation area of tidal current energy converters. Conducting a preliminary estimation of tidal current energy resources in the area is crucial for assessing the feasibility and economic viability of tidal energy power plants.

Numerical hydrodynamic models have been widely employed to assess tidal energy resources and potential environmental impacts during the early stages of tidal energy development. Gonz’alez-Gorbeña et al. (2015) performed a numerical assessment of tidal current energy potential for São Marcos Bay, Brazil, utilizing an open-source code called SisBaHiA^® [4]. Coles et al. (2017) developed a new Telemac2D model with a resolution of 1 km to assess tidal current energy at various sites around the Channel Islands [5]. This model aimed to provide a comprehensive understanding of tidal dynamics in the region and assess the potential for tidal current energy resources. The open-source hydrodynamic model ADCIRC (ADvanced CIRCulation model) was employed by Bonar et al. (2018) to assess the tidal current energy resources of multiple candidate sites in Malaysia [6]. Park et al. (2019) numerically investigated the tidal current energy resources in the Southwestern Sea of Korea using a numerical model, Modelo Hidrodinâmico (MOHID) [7]. Based on a finite-volume community ocean model (FVCOM), Karsten et al. (2008) assessed the tidal current energy in the Minas Passage, Bay of Fundy. Similarly, Wang and Yang (2020) assessed the tidal current energy in the Cook Inlet of Alaska, and Yang et al. (2020) assessed the tidal current energy in the Western Passage [8,9,10]. Chen et al. (2013) employed a 3D semi-implicit Eulerian-Lagrangian finite-element model called SELFE to explore the tidal characteristics within the Taiwan Strait and identify possible sites for the utilization of tidal current energy [11]. Based on an improvement of the originally developed SELFE model, 3D hydrodynamic model SCHISM (Semi-implicit Cross-scale Hydroscience Integrated System Model), Burić et al. (2021) investigated the potential tidal current energy resource in the strait of Novsko Zdrilo [12].

The Pearl River Estuary area has a large population, developed economy, and large energy demand [13]. Therefore, the development of tidal energy resources becomes particularly important. The Qiantang River Estuary in Zhejiang province is known to have the most abundant tidal energy resources in China, followed by the Yangtze River Estuary, Pearl River Estuary, Minjiang River Estuary, and other estuaries. Compared to the Qiantang River Estuary, the Pearl River Estuary is classified as a weak tidal estuary, and currently, there is relatively little research on tidal energy resources in this area. Therefore, the present study aims to assess tidal current energy in the Pearl River Estuary using a numerical method. The MIKE 21 Flow Model FM HD module is employed to develop a hydrodynamic model of the Pearl River Estuary, which is validated through field observation data. The model is then used to simulate the tidal currents in the Pearl River Estuary over a one-year duration. Subsequently, this study investigates the average power density and spatial distribution characteristics of tidal current energy in the region and calculates the theoretical potential of tidal energy resources in the area using the Flux method.

The paper is structured as follows: in Section 2, the numerical model for simulating the tidal currents in the Pearl River Estuary and the model validation are detailed. Section 3 introduces the tidal characteristics of the Pearl River Estuary. In Section 4, the results of the average tidal current energy power density and the theoretical potential of tidal energy resources are presented and discussed. Finally, the conclusions are presented in Section 5.

2. Materials and Methods

2.1. Model Setup

The MIKE 21 Flow Model FM (MIKE 21 FM), developed by the Danish Hydraulic Institute (DHI), is a 2D hydrodynamic model designed for simulating tidal propagation, coastal currents, wave-structure interactions, sediment transport, and water quality dynamics. MIKE 21 FM uses a flexible mesh approach, allowing elements to have variable shapes and sizes over the model domain, with triangles and quadrilateral elements. The model is based on the depth-averaged, incompressible, Reynolds-averaged Navier-Stokes equations and includes continuity, momentum, temperature, salinity, and density equations. Additionally, a Boussinesq assumption is applied, modeling turbulent eddies using an eddy viscosity. The spatial discretization of the equations employs a cell-centered volume method (DHI 2017) [14]. MIKE 21 FM has been identified as the most computationally efficient tool for accurately predicting hydrodynamic conditions in a complex estuarine environment (Symonds et al., 2016) [15]. As a result, in this study, MIKE 21 FM is used for predicting tidal propagation in the Pearl River Estuary, China.

Within MIKE21 FM, the 2D shallow water equations are obtained by integrating the horizontal momentum equations and the continuity equation over depth h = η + d. These equations are given as follows:

$$\frac{\mathsf{\partial }h}{\mathsf{\partial }t}+\frac{\mathsf{\partial }h\overline{u}}{\mathsf{\partial }x}+\frac{\mathsf{\partial }h\overline{v}}{\mathsf{\partial }y}=hS$$

(1)

$$\begin{array}{l}\frac{\mathsf{\partial }h\overline{u}}{\mathsf{\partial }t}+\frac{\mathsf{\partial }h{\overline{u}}^2}{\mathsf{\partial }x}+\frac{\mathsf{\partial }h\overline{v} \overline{u}}{\mathsf{\partial }y}=f\overline{v}h-gh\frac{\mathsf{\partial }\eta }{\mathsf{\partial }x}-\frac{h}{{\rho }_0}\frac{\mathsf{\partial }{p}_a}{\mathsf{\partial }x}-\frac{g{h}^2}{2{\rho }_0}\frac{\mathsf{\partial }\rho }{\mathsf{\partial }x}+\frac{{\tau }_{sx}}{{\rho }_0}-\frac{{\tau }_{bx}}{{\rho }_0}\\ -\frac{1}{{\rho }_0}\left(\frac{\mathsf{\partial }{s}_{xx}}{\mathsf{\partial }x}+\frac{\mathsf{\partial }{s}_{xy}}{\mathsf{\partial }y}\right)+\frac{\mathsf{\partial }}{\mathsf{\partial }x}\left(h{T}_{xx}\right)+\frac{\mathsf{\partial }}{\mathsf{\partial }x}\left(h{T}_{xy}\right)+h{u}_{s}S\end{array}$$

(2)

$$\begin{array}{l}\frac{\mathsf{\partial }h\overline{v}}{\mathsf{\partial }t}+\frac{\mathsf{\partial }h\overline{u} \overline{v}}{\mathsf{\partial }x}+\frac{\mathsf{\partial }h{\overline{v}}^2}{\mathsf{\partial }y}=-f\overline{u}h-gh\frac{\mathsf{\partial }\eta }{\mathsf{\partial }y}-\frac{h}{{\rho }_0}\frac{\mathsf{\partial }{p}_a}{\mathsf{\partial }y}-\frac{g{h}^2}{2{\rho }_0}\frac{\mathsf{\partial }\rho }{\mathsf{\partial }y}+\frac{{\tau }_{sy}}{{\rho }_0}-\frac{{\tau }_{by}}{{\rho }_0}\\ -\frac{1}{{\rho }_0}\left(\frac{\mathsf{\partial }{s}_{yx}}{\mathsf{\partial }x}+\frac{\mathsf{\partial }{s}_{yy}}{\mathsf{\partial }y}\right)+\frac{\mathsf{\partial }}{\mathsf{\partial }x}\left(h{T}_{xy}\right)+\frac{\mathsf{\partial }}{\mathsf{\partial }y}\left(h{T}_{yy}\right)+h{v}_{s}S\end{array}$$

(3)

where t is the time; x and y are Cartesian coordinates; η is the surface elevation; d is the still water depth; h = η + d is the total water depth; $\overline{u}$ and $\overline{v}$ are the velocity components along the depth-averaged and determined by $\overline{u}=\frac{1}{h}{\displaystyle {\int }_{-d}^{\eta }udz}$, $\overline{v}=\frac{1}{h}{\displaystyle {\int }_{-d}^{\eta }vdz}$; f is the Coriolis force parameter, and the expression is $f=2\mathrm{\Omega }\sin \varphi$, where Ω is the angular rate of Earth’s rotation and the value is 0.729 × 10⁻⁴ rad/s; ϕ is the geographic latitude; g is the Earth’s gravitational acceleration; ρ is the water density; ρ₀ is the reference density of water; τ_sx and τ_sy are the surface wind stress components; τ_bx and τ_by are the bottom stress components; s_xx, s_xy, s_yx, and s_yy are the radiation stress components; p_a is the local atmospheric pressure; S is the source term; u_s and v_s are the velocity components of the source term;$T_{xx}=2A\frac{\mathsf{\partial }\overline{u}}{\mathsf{\partial }x}$, $T_{xy}=T_{yx}=A\left(\frac{\mathsf{\partial }\overline{u}}{\mathsf{\partial }y}+\frac{\mathsf{\partial }\overline{v}}{\mathsf{\partial }x}\right)$ and $T_{yy}=2A\frac{\mathsf{\partial }\overline{v}}{\mathsf{\partial }y}$ are the transverse stress components; A is the horizontal eddy viscosity.

The open-boundary astronomical tidal levels are sourced from the global tidal wave model TPXO6 (Egbert and Erofeeva, 2002) [16]. Ten tidal components are utilized to compute the actual astronomical tidal processes in the deep sea. This set includes the eight main tidal components: M2, S2, K1, O1, N2, P1, K2, and Q1, and two long-period tidal currents, M_f and M_m. The hydrostatic level at the open boundary is determined based on the average sea level. The tidal level at the open boundaries is calculated by the following formula:

$${\zeta }_0(x)={\zeta }_p(x)+{\displaystyle \sum _{i=1}^{10}{A}_i}(x)\cos (w_{i}t+{\alpha }_i(x))$$

(4)

where ζ₀ represents the tidal level at the boundary, and ζ_p represents the hydrostatic level at the boundary; the index i ranges from 1 to 10, with its value corresponding to the aforementioned tidal components; $w_i$ denotes the angular frequency of the tide; A_i and α_i stand for the amplitudes and phase angles of the tide at three boundaries, respectively.

The finite volume method is used to discretize the spatial domain; that is, the space is subdivided into non-overlapping cells/elements. In the numerical simulation of two-dimensional tidal current, the shape of elements has triangles and quadrilateral elements. This meshing approach is known as unstructured grid generation technology. Its advantage is that the fitting degree of the boundary is better, and the local area can be encrypted. This model uses triangular mesh and local encryption. Two methods are used for time integration: the low-order algorithm (first-order explicit Euler method) and the high-order algorithm (second-order Runge-Kutta method). In the case of ensuring calculation accuracy, this paper employs the low-order fast calculation method.

The time step in the simulation is set to 0.01 s, with a maximum time step of 60 s. The bottom friction in the model is represented using the Manning coefficient, which is related to the bottom roughness, and its value is primarily determined by factors such as water depth, composition, and the size of the bottom sediment. Consequently, in the deep-water region of the open sea, the Manning coefficient varies within the range of 0.011 to 0.012, and in the coastal region, it typically ranges from 0.012 to 0.025 [17]. The horizontal eddy viscosity coefficient is estimated using the Smagorinsky formula, with a corresponding Smagorinsky coefficient of 0.28 m²/s.

2.2. Computation Domain and Bathymetry

The computation domain of this model includes the eight mouth gates of the Pearl River Estuary. The northern boundary of the computational domain extends from Huangpu Bridge to a point 70 m beyond Dawanshan Island. The western boundary is defined from Shangchuan Island to a longitude of 114.5° east. The domain is approximately bounded by 21.5°–41° N latitude and 116.5°–127° E longitude. The computational area has dimensions of 150 km × 160 km, with three outer sea open boundaries in the west, east, and south. An unstructured triangular grid is utilized to establish suitable boundaries and facilitate local refinement. The grid density is higher within the study area, where the grid size is approximately 1000 m, while outside the study area, the grid density is lower, with a grid size of around 2500 m. The model consists of 18,834 nodes and 29,484 elements. The model refers to mean sea level, and the underwater topographic map for the study area is based on recently published charts. The calculation domain, grid, and underwater topography are shown in Figure 1.

2.3. Model Validation

In order to validate the accuracy of the model, the numerical results are compared with the observed data. Ten tidal locations within the Pearl River Estuary were chosen to conduct measurements of tidal levels spanning from 11 April to 10 May 2011. Concurrently, the measurement of current velocities and directions was conducted at 11 stations during the same period. The spatial distribution of the measured stations is presented in Figure 2. Detailed information regarding the measured stations can be found in

Table 1. List of measured stations for tidal levels.

Station Full Name	Longitude (°, E)	Latitude (°, N)
Neilingding	113.8017	22.4265
Jinxinggang	113.6151	22.3816
Chiwan	113.8698	22.4710
Dachan’gang	113.8484	22.5452
Zhengqiang Port	113.7765	22.6570
Shanbanzhou	113.6605	22.7125
Wanqingsha	113.6287	22.5689
Hengmen	113.5199	22.5762
Nansha	113.5615	22.7435
Xianwujiao	113.6157	22.7997

and

Table 2. List of measured stations for tidal currents.

Station	Longitude (°, E)	Latitude (°, N)
S1	113.6934	22.7107
S2	113.7474	22.6201
S3	113.8019	22.5926
S4	113.7804	22.5168
S5	113.8586	22.5150
S6	113.8739	22.5333
S7	113.7227	22.4460
S8	113.8678	22.4491
S9	113.9282	22.4612
S10	113.7231	22.3392
S11	113.7952	22.3421

. Utilizing a pressure-based self-recording water level meter, water level measurements were obtained. The observation error of this measuring meter remains under 3 cm. To ensure the precision of flood and ebb tidal level recordings, meticulous hourly observations were undertaken. Furthermore, observations were executed every 10 min before and after flood and ebb tides. The measurement of tidal currents was performed utilizing an Acoustic Doppler Current Profiler (ADCP), with the profiling layer’s width adjusted between 0.25 and 1.0 m based on the water depth. Observations were executed at hourly intervals, with a minimum duration of 5 min.

The computed time series of tidal levels are compared to the observed data, illustrated in Figure 3. Similarly, the numerical results for the time series of the current are compared to the observed data, as presented in Figure 4. Negative values on the Cartesian coordinate system indicate the flood current, while positive values signify the ebb current. The model’s simulated results for tidal levels at each station correspond well with the measured results. Both high and low tidal levels are accurately predicted by the model.

To calculate the errors of current speed, the absolute relative error (ARE) was used, as shown in the following formula:

$$ARE(\%)=\left|\frac{\mathrm{Model}-\mathrm{Observation}}{\mathrm{Observation}}\right|\times 100\%$$

(5)

Error statistics of the mean current speed predictions are listed in

Table 3. Error statistics of the mean current speed between measurements and predictions.

Station	Observation (m/s)	Model (m/s)	ARE (%)
S1	0.49	0.35	−29
S2	0.33	0.41	24
S3	0.33	0.30	−9
S4	0.45	0.45	0
S5	0.41	0.52	27
S6	0.16	0.13	−18
S7	0.42	0.39	−9
S8	0.38	0.45	18
S9	0.21	0.20	−4
S10	0.45	0.39	−12
S11	0.45	0.42	−7

. From Figure 4 and Table 3, it can be observed that there are significant discrepancies between the calculated and measured values at a few measurement points. This is possibly due to differences between the actual bathymetry and the model’s bathymetry. The water depth of the region is not precisely represented at S1, S2, and S5. However, with the exception of a few measurement points, the calculated values of current speed and direction at most measurement points match well with the measured values, showing small phase deviations. The calculated values for the average speed show discrepancies of approximately 10% when compared to the measured values. Considering the overall agreement between the calculated values and the measured values for tidal levels, current speed, and direction, it can be concluded that the model provides reasonable results, accurately reflecting the tidal characteristics of the Pearl River Estuary.

3. Tidal Characteristics

3.1. The Current Speed Fields of the Rapid Flood and Ebb Tides

The tides in the Pearl River Estuary are influenced by the South China Sea tidal waves and belong to the category of irregular semi-diurnal tides. When tidal waves pass into the Lingdingyang from the outer sea, the flared bay shape causes the accumulation of tidal energy along the way, resulting in an increasing tidal difference from the mouth of the bay to the top of the bay.

Based on the results of numerical simulations conducted on a large scale, Figure 5 presents the velocity fields of the maximum flood and ebb tidal currents during spring tides. From the figure, it can be observed that the current speed in the main channel is stronger than that in the shoals, with the eastern area having higher velocities than the western area. The mainstream of the flood and ebb currents tend to align with the main channel, and there is a noticeable diverging current at the intersection of the waterways between the west shoal and the main deep channel. At the moment of the rapid ebb tide, the tide flows from the upstream waterways and discharges into Lingdingyang together. The falling tide in the east trough of Lingdingyang moves southward and continues its southward flow after merging with the falling tide in Shenzhen Bay. In the Tonggu Sea, it splits into two strands at the top of Hong Kong Airport. One strand flows eastward out of Lingdingyang through the Hong Kong waterway, while the other turns southwestward and flows southeastward around Lantau Island. At the moment of rapid flood tide, the tide flows from the east to the west, entering the Lingdingyang Estuary through the Hong Kong Channel. It combines with the rising tide from the Zhuhai to Lantau Island section in the Tonggu Sea and then it moves northward. In the vicinity of Chiwan, the current bifurcates, with one branch entering Shenzhen Bay and the other persisting in its northward trajectory.

In general, the Pearl River Estuary is a weak tidal estuary characterized by small tidal ranges and low current velocities. On average, the current speed during the flood tide is about 0.5 m/s, while during the ebb tide, it is around 0.3 m/s.

3.2. The Annual Current Speed

Figure 6 and Figure 7 display the distribution of annual average current speed and annual maximum current speed, respectively. To comprehensively understand the theoretical potential of tidal current energy resources in the Pearl River Estuary, three cross sections have been selected based on the available bathymetric data. From the figures, it can be observed that the current speed in the main channel is stronger than that in the shoals, and the current velocities are generally higher in the eastern area compared to the western area, which is similar to the characteristics of flood and ebb tidal currents. The annual average current speed for Cross section I, Cross section II, and Cross section III are 0.35 m/s, 0.37 m/s, and 0.33 m/s, respectively. The annual maximum current speed for Cross section I, Cross section II, and Cross section III are 0.94 m/s, 0.85 m/s, and 0.65 m/s, respectively. The annual average current speed in the main channel is approximately 0.5 m/s, while the annual maximum current speed is approximately 1.2 m/s. These results indicate a weak tidal characteristic in the study area.

4. Tidal Current Energy Estimation

4.1. Average Tidal Current Energy Power Density

One objective of this paper is to assess the potential of the tidal current resource in the Peal River Estuary in terms of power density. The power density per square meter is proportional to the cube of the current speed.

$$P=\frac{1}{2}\rho V^3$$

(6)

where P is the tidal power density (W/m²), ρ is the seawater density (kg/m³), and V is the tidal current speed (m/s).

The monthly average power densities in the Pearl River Estuary for January (winter), April (spring), July (summer), and October (autumn) are shown in Figure 8. Consistent with the current speed distribution, the typical monthly tidal energy power density in the main channel is higher than that on the shoals. For Cross section I, the monthly average power densities in January, April, July, and October are 0.05 kW/m², 0.05 kW/m², 0.05 kW/m², and 0.05 kW/m², correspondingly. In Cross section II, the average monthly power densities for January, April, July, and October are 0.06 kW/m², 0.05 kW/m², 0.06 kW/m², and 0.06 kW/m², respectively. Within Cross section III, the monthly average power densities for January, April, July, and October amount to 0.04 kW/m², 0.04 kW/m², 0.04 kW/m², and 0.04 kW/m², respectively. The average power densities of tidal currents in summer and autumn are larger, followed by winter. The spring season has the lowest average power density.

The numerical calculations were recorded every half an hour, resulting in the generation of flow fields for each time interval. Based on these flow fields, the corresponding tidal current power density was calculated for each specific time. Finally, by averaging these values, the distribution of annual average tidal current power density was obtained in Figure 9.

Figure 9 reveals that the distribution of annual average tidal current power density in the Pearl River Estuary is generally consistent with the distribution of the tidal current fields. The overall average tidal current power density is relatively small, with values mostly below 0.10 kW/m².

4.2. Potential Resources

Based on the Flux method (Lu and Qiao, 2008) [18], the theoretical potential of tidal current energy resources, P_total, is calculated using the following formula:

$$P_{total}=P_{\mathrm{m}}\cdot A$$

(7)

where P_total represents the theoretical average power of tidal current energy in kilowatts (kW), P_m represents the average power density of tidal current energy in kilowatts per square meter (kW/m²), and A represents the cross-sectional area of the waterway in square meters (m²).

The characteristics of tidal current power density distribution are shown in

Table 4. Statistics for reserves in the theory of tidal current energy in the estuary of Pearl River.

Cross Section	Width (m)	Average Water Depth (m)	Cross-Sectional Area (m2)	Annual Average Power Density (kW/m2)	Theoretical Potential Resource (10,000 kW)
I	26,466	8.4	223,110.7	0.05	1.1
II	28,455	6.4	183,252.5	0.05	0.9
III	29,187	4.5	129,883.7	0.04	0.5

. For Cross section I, the annual average power density of tidal current energy is 0.05 kW/m², resulting in a theoretical potential of 11,000 kW. For Cross section II, the annual average power density of tidal current energy is 0.05 kW/m², resulting in a theoretical potential of 9000 kW. For Cross section III, the annual average power density of tidal current energy is 0.04 kW/m², resulting in a theoretical potential of 5000 kW. Section I has the highest theoretical potential resource, followed by Section II, and Section III has the lowest potential. The theoretical potential resources tend to increase with the increase in water depth.

5. Conclusions

In the present work, an unstructured-grid MIKE 21 FM tidal hydrodynamic model was built to simulate tidal hydrodynamics and assess the tidal energy potential in the Pearl River Estuary. The numerical model was validated by using the field observation data. The model effectively reproduces tidal hydrodynamics within the study area. The summarized conclusions are as follows:

• The distribution of annual average tidal current power density in the Pearl River Estuary is generally consistent with the distribution of the tidal current fields. The average power densities of tidal currents in summer and autumn are larger, followed by winter, and it is the smallest in spring.
• The annual average power density of tidal energy is generally smaller than 0.10 kW/m². The theoretical resource potential increases with the increase in water depth. The theoretical potential of tidal energy resources in the Pearl River Estuary was finally assessed to be about 11,000 kW.
• The tidal range in the Pearl River Estuary is small, resulting in a relatively weak tidal power and low average power density of tidal energy. Therefore, the tidal energy resources in the estuary are limited.