Keywords

1 Introduction

Coastal zones have always been preferential areas for human settlement mainly because they offer essential resources, but also spaces for recreational and cultural activities. They offer points of access to world-wide marine trade and transport. Consequently, population density is significantly higher in coastal areas [1] and the migration to coastal zones is expected to continue into the future [2]. Population rise in these zones also increases the number of people exposed to the risks that characterize these areas, such as extreme natural events [3] and floods due to sea-level rise [4].

Human occupation and development in the coastal zones generates high pressure on the natural resources and ecosystems [5]. A high number of activates, ranging from recreation to commercial, compete for the limited resources of the ocean and often lead to negative impacts [6]. These growing pressures require management approaches that promote a conscious use of the ocean resources, fulfilling current needs without compromising future generations and reducing impacts on the ecosystems. Due to higher spatial and temporal variability, the management of ocean and coastal areas is more demanding than terrestrial ecosystems management [7], and dynamic management approaches are necessary to match the dynamism of the system [8].

An effective dynamic management of marine systems requires real-time sources of information [9]. Numerical models are important tools and sources of information to assist in marine activities by continuously providing estimates and forecasts on the state of the ocean. In addition, models results can also be combined with other data and increase our knowledge and managing capacity of the marine space [10, 11].

This work describes the implementation of a three-dimensional (3D) hydrodynamic computational model for the coastal zone in Sines, Portugal. Results from the model will be used as a management tool for the port activities, including the aquaculture. This model implementation is part of the ongoing Project PiscisMod (16-02-01-FMP-0049), financed by Mar2020.

2 Case Study

Sines is a city located on the west coast of Portugal, with approximately 15 000 inhabitants [12] and a floating population of about 5 000 people, due to economic and touristic reasons. Sines plays a major role in terms of energy production and storage. There are two large production centers of oil and gas, GALP refinery and the Repsol petro-chemical industrial complex, both connected via pipelines to oil-bearing and petro-chemical terminal of port of Sines [13].

The port of Sines, located south from the city, is one of the most important deep-water ports of the country, having geo-physical conditions to accept a wide variety of ships. It is the country’s leading energetic supplier (crude and its derivatives, coal and natural gas) and an important gateway for containerized cargo. The port consists of five terminals – liquid bulk, liquid natural gas, petrochemical, container, and multipurpose – as well as fishing and leisure ports (Fig. 1) [14]. In recent studies the port of Sines was found to be the most efficient port in Portugal in terms of emissions per cargo ship, although it is also the most pollutant due only to the fact that it receives more cargo than other ports [15]. Near the container terminal, protected by the breakwater, there is a delimited area allocated for the development of aquaculture (Fig. 1) consisting of 16 cages for a yearly production of 500 metric tons of European sea bass (Dicentrarchus labrax). The Sines thermoelectric power plant is located around 3.5 km southeast of the aquaculture zone and uses seawater to cool the generators. The seawater intake and restitution points are protected by jetties and located close to São Torpes beach. The water is discharged by two open channels with 4.5 m depth distant about 400 m from the intake. On a yearly average the power plant system uses 40 m\(^{3}\)/s of water [16].

This coastal area, as well as the rest of the west coast of Portugal, is frequently under the influence of upwelling, caused by the predominant north wind regime [17].

Fig. 1.
figure 1

Port of Sines terminals (left panel) and a close-up on the aquaculture facilities area (right panel) with the monitoring stations (S1 to S4).

Full size image

3 Methodology

3.1 Model Description

The numerical model used in this work is the MOHID water modelling system (www.mohid.com), a model developed at the Marine and Environmental Technology Research Center (MARETEC). MOHID simulates physical and biogeochemical processes in the water column and the interactions between the water column with the sediment and the atmosphere. This model has already shown its ability to simulate complex systems and processes including those present in coastal and estuarine system [18, 19]. And it was already proven useful as a management tool in aquaculture activities [20, 21]. The model solves the three-dimensional primitive equations for the surface elevation and 3D velocity field for incompressible flows. The hydrostatic, Boussinesq and Reynolds approximation are assumed. It solves the advection-diffusion of temperature and salinity. The density is calculated using the UNESCO state equation as a function of salinity, temperature and pressure [22]. The viscosity is separated into its horizontal and vertical components, being the horizontal usually set to a constant. The vertical turbulent viscosity is handled using the General Ocean Turbulence Model (GOTM), which is integrated into the MOHID code, having its k-\(\varepsilon \) model parameterized according to Canuto et al. [23].

In the MOHID model the discretization of the equations is done using a finite volume approach. The discrete form of the governing equations is applied to a cell control volume, making the way of solving the equations independent of cell geometry and allowing the use of generic vertical coordinates [24]. Horizontally the equations are discretized in an Arakawa-C staggered grid [25]. For time discretization the model solves a semi-implicit ADI (Alternating Direction Implicit) algorithm with two time levels per iteration as proposed by Leendertse [26] and implemented as explained in Martins et al. [24]. The numerical scheme chosen for the advection and diffusion is the Total Variation Diminishing (TVD) method which consists in a hybrid scheme between a first order and a third order upwind method using a ponderation factor calculated with the SuperBee method [27].

3.2 Model Setup

The Sines modelling system presented in this study was implemented following a downscaling methodology, which is useful to provide boundary conditions from regional to local scale models [28]. A system of five nested domains of increasing horizontal resolution was configured. The first domain is a data acquisition window from the PCOMS (Portuguese Coast Operational Modeling System) model, a regional model for the Portuguese coast with a 0.06\(^\circ \) horizontal resolution and a vertical discretization of 50 layers [29]. From domain 2 to domain 5 the horizontal resolution increases by a factor of 4 (0.015\(^\circ \), 0.00375\(^\circ \), 0.000938\(^\circ \) and 0.000234\(^\circ \)), the most refined nested model (domain 5) has a horizontal resolution of approximately 25 m. The bathymetric information required for the model was obtained from the European Marine Observation and Data Network (EMODnet) [30] for domains 2 and 3 while for domains 4 and 5 the information was provided by the Instituto Hidrográfico (www.hidrografico.pt) of the Portuguese navy. The result of the interpolations of the bathymetry into the different domains grids are presented in Fig. 2.

Fig. 2.
figure 2

Model domains: location and bathymetry. Black marker represents the buoy location.

Full size image

All domains are set in 3D and follow a generic vertical discretization approach [31] allowing the implementation of 7 sigma-type layers close to the surface and below a variable number of cartesian-type layers depending on the maximum depth of the corresponding domain.

The first domain model is a data acquisition window from the PCOMS model with a high temporal resolution, 900 s, enough to represent the main processes coming from the open ocean including the tide signal. PCOMS results for water level, velocity fields, temperature and salinity will be used as boundary conditions for domain 2 using a Flow Relaxation Scheme (FRS) for the velocities, temperature and salinity [32] while the water level is radiated using the Flather method [33]. On the open boundary with the atmosphere the model computes momentum and heat fluxes that interact with the free surface using results from the Mesoscale Meteorological Model 5 (MM5) atmospheric forecast model for the west Iberian Coast, running at IST (http://meteo.tecnico.ulisboa.pt). This model provides hourly results for air temperature, wind intensity and direction, atmospheric pressure, relative humidity, solar radiation and downward long wave radiation with a 9 km spatial resolution. For domain 1 to domain 4 the simulations are performed online (all at the same time) and domain 4 provides a window of results with a 900 s interval that will be used as boundary condition for domain 5 simulations. This configuration allows running domain 5 in a different computer to improve the computational time of the system. Time-step values used for the model calculations were 60.0, 30.0 and 7.5 s for domains 2, 3 and 4 respectively and 4.0 s for domain 5. Average duration for running this setup to obtain simulation results for 1 day, on a Windows computer with an Intel© CoreTM i7-980 processor, was 1 h and 20 min from domain 1 to domain 4 and 50 min for domain 5. Summing to 2 h 10 min for the whole Sines modeling system.

The thermoelectric power plant intake and discharge are modeled in domains 3 and 4. In both domains the water intake was modeled by a sink term in the horizontal grid cell closest to the location of the real intake, a flow of 40 m\(^{3}\)/s was considered according to data from the power plant environmental declaration [16]. The discharge is modeled by a source term: in domain 4, two source terms where imposed in different horizontal grid cells in order to represent the two discharge open channels and a flow of 20 m\(^{3}\)/s was considered for each source term; in domain 3 only one source term with 40 m\(^{3}\)/s flow was considered due to the model spatial resolution. The water intake and discharges are considered along the whole water column to simulate the real conditions. A bypass function was implemented to model the 10 \(^\circ \)C rise in temperature of the water in the discharge when compared to the intake caused by the cooling of the power plant, all other properties have the same value on both sides. The temperature increase was defined according to previous studies by Salgueiro et al. [35]. The breakwater that protects the container and multipurpose port terminals (see Figs. 1 and 2) was modeled only in domain 5. The breakwater modeling was done by computing the deceleration that the obstacle causes on the fluid that passes through it by using a drag equation:

$$\begin{aligned} \frac{F}{\rho Vol} = \frac{1}{2} v^2 \,C_d \frac{A}{Vol} \end{aligned}$$
(1)

F is the drag force caused by the obstacle, \(\rho \) is the density of the water, Vol is the volume of the computational cell, A is the surface area and \(C_d\) is the drag coefficient of the obstacle. The modelling of the breakwater was done by attributing a drag coefficient of 2 to the computational cells chosen to represent the structure.

4 Results and Discussion

4.1 Model Validation

Two sets of data were used to validate the Sines modelling system. Domain 2 and 3 model results were validated with temperature and water level. Temperature observations were recorded by a moored buoy located near the Sines coastal area (37\(^\circ \)55.3’ N, 008\(^\circ \)55.7’ W, black symbol in Fig. 2) and water level observation from a tide Gauge located inside the port of Sines [36]. Model results obtained in domain 5 simulations were validated with in situ data recorded on 29/06/2018, between 12:15 h and 14:15 h, on four stations (S1, S2, S3, S4 in Fig. 1) inside the container terminal near the aquaculture area. S1 is the station most protected by the breakwater and corresponds to a boat access platform, S2, S3 and S4 are progressively closer to the end of the breakwater. S2 is a support platform for the aquaculture activities and S3 and S4 are buoys that mark the perimeter of the aquaculture area. Depth profiles for velocity, temperature, salinity and density were recorded on all stations. Velocity was recorded using an Aanderaa acoustical current meter (model DCS4100) and the water properties were recorded with a CTD from Falmouth Scientific, Inc (model FSI NXIC CTD).

Domain 2 and 3 Validation. Model validation was made comparing time series of temperature and water level. For temperature a sample size of 1998 data points obtained between 26/12/2015 and 10/07/2016 were used. The analysis of temperature results for domain 2 and 3 show a significantly small difference from the observations (Fig. 3), having a Root Mean Square Error (RMSE) of 0.62 and 0.66 \(^\circ \)C, respectively. Pearson correlation assumes values of 0.90 for domain 2 and 0.89 for domain 3 which indicates a good fit between models and observations. The Bias Error (BE) of \(-0.18\) and \(-0.21\,^{\circ }\)C for domains 2 and 3 respectively indicate that the model tends to underestimate the temperature.

For water level a sample size of 2647 data points between 01/01/2016 and 12/05/2016 were used. Model results show a significantly small difference from observed values (Fig. 3) which is shown by the RMSE values for domain 2 and 3 of 0.07 m. Pearson correlation assumes a value of 0.997 for both domains also showing a good fit between model and observations.

Fig. 3.
figure 3

Model results for domain 3 and buoy records for temperature and water level.

Full size image

Domain 5 Validation. In domain 5, first the hypothesis of the breakwater being permeable to water was tested. A scenario with the breakwater modeled with a drag coefficient of 2 (model A in Figs. 4, 5, 6 and 7) was compared with a scenario in which the breakwater is impermeable by setting the horizontal grid cells that represent it to land cells (model B in Fig. 4, 5, 6 and 7). Model results from each model were compared with in situ data from the stations inside the container port. Differences in velocity profiles are easily noted, model A shows results closer to observations in 3 of the 5 situations compared, particularly in station S1 at 11:45 h and S2 at 13:45 h and 14:00 h (see Fig. 4). On S2 at 12:45 h although observation values are closer to model B results than model A they are significantly more distanced from either model than on other occasions, this may be explained by these two data points being closer to the surface then most. Inside the port there is movement from container ships and container ship towboats, these movements might cause waves that increase the water velocity in a way that the model can’t reproduce. On S2 at 13:15 h model B is closer to observations but model A shows the same tendency of decreasing velocity with water column depth in the locations of the data points while model B shows the opposite tendency. The difference between the models for salinity were negligible (see Fig. 6). For density, model B shows, in general, lower values than model A but, on average, the deviation is not significantly higher than in model A. Notable differences between model A and B are seen for temperature with model A having results significantly closer to observation values. Model A was chosen to represent this system due to being closer to measured values.

Model A validation was done by comparing model results with in situ observations for velocity, temperature, salinity and density. RMSE, BE and Pearson correlation were calculated for temperature, salinity and density but not for velocity due to the lack of data points. For the velocity comparison, the biggest difference is found on station S2 at 12:45 h with an error of 5.39 cm/s at 2 m of depth and 4.31 cm/s at 4 m depth. Excluding these 2 data points, which as previously explained might be influenced by boat movement in the area, the next biggest error is 1.64 cm/s at 14.5 m depth on S2 at 13:45 h. Overall model results are close to observations with a tendency to underestimating reality. Regarding the comparison between temperature model results and observed data for the stations S1, S2, S3 and S4 the RMSE values are 0.33, 0.65, 0.58 and 0.75 \(^\circ \)C respectively, the BE are 0.22, 0.65, 0.55 and 0.69 \(^\circ \)C and the Pearson correlations are 0.96, 0.96, 0.77 and 0.98. The difference between model results and observations are relatively small with a tendency for overestimation by the model, correlation shows a good fit but on location S3 this is significantly worse than on other stations due to the different temperature gradient of the observations when compared with the behavior of the gradient on other stations (see Fig. 5). For salinity, the RMSE values range from 0.14 to 0.16 PSU – 78 and the BE assumes the same value as the RMSE for each location, this indicates that the model overestimates salinity on all compared data points, which is confirmed by observing the depth profiles in Fig. 6. Correlation values for stations S1 to S4 are \(-0.57, 0.66, 0.73\) and 0.72, respectively, indicating a moderate fit for the data except on station S1, this value can be explained by a low number data points in S1, 12 data points, which might not be enough to accurately represent the trend that salinity follows with increasing depth, all other observations depth profile have more than 20 data points and show an increase of salinity with depth while S1 shows a decrease. Density RMSE values are 0.12, 0.04, 0.09 and 0.06 kg/m\(^{3}\), BE values are 0.11, 0.03, 0.07 and 0.01 kg/m\(^{3}\) and correlation values are 0.97, 0.96, 0.82 and 0.98, all from S1 to S4 respectively. Statistical indicators show a good simulation of the density by the model with the biggest difference being shown in S1. This was expected due to density results being directly dependent on temperature and salinity results. As these water properties were relatively accurately modeled, density should be too. Density depth profiles are presented in Fig. 7.

Fig. 4.
figure 4

Depth profiles for velocity: model results and observations.

Full size image
Fig. 5.
figure 5

Depth profiles for temperature: model results and observations.

Full size image
Fig. 6.
figure 6

Depth profiles for salinity: model results and observations.

Full size image
Fig. 7.
figure 7

Depth profiles for density: model results and observations.

Full size image

Looking at model results for surface velocity for domain 4 (Fig. 8) it is possible to confirm that the breakwater is having the desired effect of protecting the inside of the port from strong currents, while still allowing some passage of water. The strongest currents occurred from 14:00 h to 17:00 h, coincident with strong winds in the same direction.

Fig. 8.
figure 8

Model results for velocity (intensity and direction) on 29/06/2018.

Full size image

The deviations between model results and observations may come from a variety of sources. The model relies on results from previous models, other hydrodynamics models and atmospheric models; errors from these models will propagate to the model in this study. While the variation of properties in the water is continuous in space and time in the model, values are represented discretely causing a loss of information specially when high variability exists. The breakwater is, in fact, larger in width at the bottom of the water than at the surface but due to lack of information it was modeled as having the same width and the chosen drag coefficient might not be adequate. As such, more validation data is required to assess it.

5 Concluding Remarks

Model results show a good representation of the hydrodynamics and water properties in the study area. Knowing this, model results can prove useful for management of activities inside the container port, specially the aquaculture. Future work from this project will focus on this. In the next stage of this project the model will be compared with more data for further validation for different conditions. To improve accuracy, the weather forecasting model currently in use (MM5) will be replaced by a WRF (Weather Research and Forecasting) model for the Portuguese area. This model offers a better spatial resolution of 3 km, when compared to the 9 km from the MM5 model.

In the future, more processes will be added to model the water quality of the area and a bio-energetic module using DEB (Dynamic Energy Budget) theory [37] will also be included. This will allow the modelling of the fish aquaculture in order to optimize the feeding process and to asses the impact of this activity on the local ecosystem.