Assessing the Impact of Chemical Loads from Agriculture Holdings on the Puck Bay Environment with the High-Resolution Ecosystem Model of the Puck Bay, Southern Baltic Sea

Abstract

This paper describes the ecohydrodynamic predictive model EcoPuckBay—the ecosystem part—for assessing the state of the Puck Bay coastal environment and its ecosystem. We coupled the EcoPuckBay model with the land water flow models (Soil and Water Assessment Tool (SWAT) for surface water and Modflow for groundwater). To evaluate the quality of the results obtained from the EcoPuckBay model, a set of basic statistical measures for dissolved oxygen, chlorophyll-a, nitrates, and phosphates were calculated, such as mean, Pearson correlation coefficient (r), root-mean-square-error (RMSE), and standard deviation (STD). The analysis presented in this paper shows that the EcoPuckBay model produces reliable results. In addition, we developed a nutrient spread module to show the impact of agricultural activity on the waters of the Puck Bay. The EcoPuckBay model is also available in operational mode where users can access 60-h forecasts via the website of the WaterPUCK Project through the “Products” tab.

1. Introduction

Puck Bay is an example of a region that is highly vulnerable to anthropogenic impact. Therefore, it has been included in Natura 2000 (https://ec.europa.eu/environment/nature/natura2000). As a result, it requires preservation or restoration of the favourable conservation status of species and habitats by introducing appropriate protection measures. The strategic actions and the policy of the authority of the Puck District regarding the environmental protection involve not only the respect of the Natura 2000 legislation but also the execution of European legislation including Water Framework Directive, Marine Strategy Framework Directive, Habitats Directive, Baltic Sea Action Plan, and the strategic program of the environment protection for the Puck Commune. The main aims of the Puck Commune policy are the improvement of the environment, sustainable development, resilience to climate change, and protection of natural resources such as water.

WaterPUCK service, developed as a part of the WaterPUCK project (https://www.waterpuck.pl/en/), is a tool to carry out analyses and forecasts and to support key decisions for the management of areas where there is or hypothetically could be agricultural activity. WaterPUCK allows determining the impact of these activities on the surface and underground waters in the area of the analysed municipality and on the coastal waters of the Puck Bay.

The EcoPuckBay model has been developed as part of the WaterPUCK project. The hydrodynamic part of the EcoPuckBay has been described and validated [1]. The aim of the project is to create integrated information and predictive service for the Puck Commune through the development of a computer system providing the WaterPUCK service, which will clearly and practically assess the impact of farms and land use structures on surface waters and groundwater in the Puck Commune and consequently on the quality of the waters of the Puck Bay [2,3]. The construction of the service is based on in situ research, surveys, environmental data (chemical, physicochemical, and hydrological), and numerical modelling. WaterPUCK service is an integrated system consisting of computer models interconnected with each other, operating continuously by supplying it with meteorological data and combines four main modules:

• two calculators for farms in the Puck Commune as an interactive applications [4,5],
• a numerical model of groundwater flow based on Modflow [6],
• a comprehensive model of surface water runoff based on Soil & Water Assessment Tool (SWAT) [7],
• a three-dimensional numerical model of the Puck Bay ecosystem consisting of a hydrodynamic [1] and biochemical part with a nutrient spread module (this paper).

The spatial distributions of major biochemical tracers (e.g., nutrients, dissolved inorganic carbon, and oxygen) are governed by a combination of physics, chemistry, and biology. High-resolution numerical ecosystem models play an important role in deciphering coastal nutrient dynamics, in providing a quantitative framework for assessing the contributions of different processes, and in interpreting field observations. The following sections of the paper present the study area, the configuration of the EcoPuckBay model, its linkage with the surface water flow models and groundwater flow model, statistical analysis of the obtained results compared with both environmental and numerical data, an introduction of the nutrient spread module, and the operational mode of the model with 60-h forecasts.

2. Materials and Methods

2.1. Study Area

The southern part of the Baltic Sea enclosing the Puck Commune is a popular tourist region that is also heavily influenced by an anthropogenic activity of local residents and farming. This makes the Puck Bay a natural reservoir for waste deposition of fertilisers and other inputs delivered through the soil, groundwater, river, or direct deposition.

This region is massively shaped and influenced by several factors that can pressure the state of physical and biochemical parameters. One of the most important factors influencing the region’s unique ecosystem is topography. The average depth of the Gulf of Gdańsk is about 50 m, with a maximum depth (Gdańsk Deep) of 118 m. From the northeast, it is surrounded by Hel Peninsula, which serves as a natural barrier for mixing with the open waters of the Baltic Sea, keeping the salinity ranging mostly within the range of 7–8 PSU with a deviation of around 1 PSU. Puck Bay is also heavily influenced by the river discharge from land, resulting in lowered salinity, especially in the coastal surface waters. The largest river in the region is the Vistula (for which the mouth is denoted as a symbol “1” in the left-hand side of Figure 1), discharging an average of over 1000 m${}^3$ s${}^{-1}$ of freshwater. Water temperature in the region ranges from over 20 ${}^{\circ }$C at the surface during summer, with the maximum usually in August, to around 2 ${}^{\circ }$C in February. Water stratification is frequent during the warmer months, leading to occurrence of seasonal thermoclines. During the winter seasons, thermocline declines and the water becomes well mixed.

In order to assess the possibility and scale of eutrophication and water pollution, the area of interest and effective domain (for which size is reduced compared to the previous version [1] to focus mainly on changes occurring closer to the Puck Commune, where the strong influence of the waters flowing with the Vistula River is no longer as great as near its mouth) in the EcoPuckBay model covers the western part of Gulf of Gdańsk (Figure 1). It can be divided further into a shallow part known as Puck Bay and the semi-enclosed Puck Lagoon to the northwest.

2.2. EcoPuckBay Configuration

The EcoPuckBay model originates from Community Earth System Model (CESM)-coupled global climate model by National Center for Atmospheric Research (http://www.cesm.ucar.edu/models/ccsm4.0). CESM is a state-of-the-art model system consisting of five separate components with an additional coupler controlling time, exciting forces, domains, grids, and information exchange between the models. For the purpose of the WaterPUCK project, CESM was downscaled and adapted for the Puck Bay region for further development at the Institute of Oceanology, Polish Academy of Sciences.

The EcoPuckBay ecosystem module is based on the nutrient–phytoplankton–zooplankton–detritus (NPZD) approach [8]. The model predicts nutrient distributions (N, Si, and P), three phytoplankton functional types (diatoms, picophytoplankton/nanophytoplankton, and diazotrophs (nitrogen-fixing organisms)), chlorophyll-a as separate variables, zooplankton, pelagic detritus, dissolved oxygen, and pesticides (glyphosate, diflufenican, metazachlor, chlorpyrifos, and antrachinon). Sources of nutrients include atmospheric deposition and sedimentary sources. Sinking particles can be associated with ballast, and the presence of ballast changes the remineralization length scale assumed for the model. Many, though not all, models can be cast in a general form as a coupled set of time-dependent advection, diffusion, reaction equations:

$$\frac{\partial S}{\partial t}+(V+w_s)\nabla S-\sum _{i=1}^3\frac{\partial }{\partial x_i}\left(K_{x_i}\frac{\partial S}{\partial x_i}\right)=F_S,$$

(1)

where S refers to a set of prognostic or predicted biochemical variables. The second and third terms on the left side of the equation describe the advection and mixing, respectively, where $V(u,v,w)$ is the velocity vector, $w_s$ is the sinking velocity of pelagic detritus, and $K_{x_i}$ is a turbulent diffusion coefficient. All of the chemical and biological interactions are depicted as one term, $F_S$, on the right-hand side, and the specific equations are listed in the 3D CEMBS model description [9].

2.3. Land–Water Linkage

We coupled the EcoPuckBay model from the land side with two models SWAT (surface water) and Modflow (underground), which we briefly describe below.

2.3.1. Surface Water

Information about the water volume discharged by rivers with their mouths located within the Puck Commune (Figure 1 and

Table 1. Locations of mouths of watercourses included in the EcoPuckBay model domain.

	Watercourse	Longitude	Latitude
1	Vistula	18.95	54.35
2	Bold Vistula	18.78	54.37
3	Still Vistula	18.66	54.41
4	Oliwski Stream	18.60	54.42
5	Kamienny Stream	18.56	54.46
6	Kacza	18.56	54.48
7	Ściekowy Canal	18.51	54.61
8	Reda	18.47	54.64
9	Gizdepka	18.46	54.66
10	Potok Błądzikowski	18.45	54.70
11	Płutnica	18.39	54.72

, numbers from 8 to 11) is being provided by the hydrological model SWAT that has been implemented as one of the WaterPUCK project’s stages [7]. The SWAT model includes the preparation of the innovative and complex hydrological model coupled with the nutrient concentration module including meteorological data (precipitation, wind, temperature, and atmospheric pressure). The proposed solution is based upon real-time observation (local weather station) and on short-term weather forecasts (the Interdisciplinary Centre for Mathematical and Computational Modelling University of Warsaw (ICM-UW) website). The hydrological computations are performed with SWAT software. The transformation of precipitation data into surface runoff have been achieved with the SCS (Soil Conservation Service) curve number procedure through the accumulated runoff volume and the time of concentration (the time from the beginning of a rainfall event until the entire subbasin area contributes to flow at the outlet) [7].

Information on the deposition of nutrients from the HYPE model is available in the form of monthly averages. As a result of the HELCOM directives in force, the actual values coming into the Baltic Sea from Polish areas have been strongly reduced over the last 30 years [10]. Therefore, the use of 30-year averages would lead to an overestimation and distortion of the actual runoff. Therefore, the current deposits for the HYPE rivers (numbers 1 to 7, Figure 1) were established on the basis of [10]. On the basis of the abovementioned work, actual concentrations were determined for nitrates (0.9 mg dm${}^{-3}$), ammonia (0.07 mg dm${}^{-3}$), phosphates (0.07 mg dm${}^{-3}$), and silicates contained in the waters of river origin entering the Baltic Sea (1.1 mg dm${}^{-3}$). The concentrations thus established were related to the average daily volumes of freshwater introduced by these rivers, resulting in satisfactory deposition rates close to real ones.

The SWAT model operates in operational mode, producing forcing data for the EcoPuckBay model for the land/water border; however, its area of operation does not cover the entire domain of the EcoPuckBay effective model. SWAT provides information only on the rivers within the Puck Commune, while the EcoPuckBay model also includes the rivers south of the Reda River, i.e., Ściekowy Canal, Kacza, Kamienny Stream, Oliwski Stream, Still Vistula, Bold Vistula, and the Vistula. Therefore, for full information on freshwater streams in the EcoPuckBay model, this data had to be obtained from another existing hydrological model. For this purpose, numerical values from HYPE model simulations based on historical time series (available: 1980–2010) are used. Its geographical domain covers the whole of Europe, and public data is available in a time resolution analogous to the SWAT model, i.e., as average daily flows in m${}^3$ s${}^{-1}$. Four subcatchments were separated from the HYPE model—one for each of the Vistula mouths and one for the remaining rivers on the basis of a uniform distribution.

Among the river data obtained from the HYPE model, the Vistula is the most important source of both the volume of freshwater flowing into the Gulf of Gdansk and the deposition of nutrients (

Table 2. The river runoff into the Baltic Sea in 2016 inside the model domain border.

	Watercourse	Mean Runoff	Nitrogen	Ammonia	Phosphorous
		[m${}^3$ s${}^{-1}$]	[tons year${}^{-1}$]	[tons year${}^{-1}$]	[tons year${}^{-1}$]
1	Vistula	1064.20	15,102.32	1174.63	1174.63
2	Bold Vistula	2.06	29.17	2.27	2.27
3	Still Vistula	6.07	86.08	6.69	6.69
4	Oliwski Stream	0.35	4.98	0.39	0.39
5	Kamienny Stream	0.35	4.98	0.39	0.39
6	Kacza	0.35	4.98	0.39	0.39
7	Ściekowy Canal	0.35	4.98	0.39	0.39
8	Reda	0.28	36.33	3.24	0.94
9	Gizdepka	0.49	29.74	1.83	0.72
10	Potok Błądzikowski	0.31	25.17	0.83	0.36
11	Płutnica	1.35	32.62	9.04	3.79

), which has a huge impact on environmental conditions within the effective domain of the model. The order of magnitude of the Vistula’s water flow is on average about three orders of magnitude larger than the other freshwater sources in the EcoPuckBay model. The annual cycle of the Vistula has been similar in character for decades, with peaks in the spring months—which is related to snow melt and summer minimum. Thanks to this, the use of historical time series for the Vistula is a good approximation of the actual flow for the last few years.

2.3.2. Groundwater

A numerical transport model based on the Modflow code developed at another stage of the project allows determining the nitrate load in the groundwater flowing into the Puck Bay. The MT3DMS numerical code was applied to solve the advection–reaction–dispersion equation [11]. Transient calculations were performed with the third-order total-variational diminishing (TVD) numerical method for the advective term, while for the steady-state solution, the standard finite difference discretization was applied. This model was calibrated based on actual NO${}_3$ concentration values and was joined with the EcoPuckBay model through a coupling module.

2.4. Water–Water Border

The results of the 3D EcoPuckBay model are limited to the effective area of Puck Bay. However, the entire model grid covers a wider area marked in Figure 1. This is to ensure that boundary conditions are properly simulated. Along the line of the northern border of the 3D EcoPuckBay model, data from the 2.3 km 3D CEMBS prediction model are transferred to the EcoPuckBay model (see Figure 1). Results from 3D CEMBS [9,12] are used to provide forcing fields in the EcoPuckBay model through sequential information transfer. The mechanism of this module is to interpolate values from 3D CEMBS to EcoPuckBay model’s grids.

2.5. Data Used for Evaluation

To evaluate EcoPuckBay model, we used in situ measurements taken throughout the monitoring activities of the Voivodship Inspector of Environmental Protection in Gdańsk as well as the numerical data of biochemical parameters calculated with the ice-ocean model NEMO-Nordic (Nucleus for European Modelling of the Ocean) coupled with the biogeochemical model SCOBI (Swedish Coastal and Ocean Biogeochemical model) acquired from the Marine Copernicus database. A brief summary description of the data used for the evaluation in this paper is presented below.

2.5.1. VIEP

The biggest set of the in situ data used for the EcoPuckBay model evaluation performed in this paper is samples collected during regular measurements of the marine environment of the Baltic Sea conducted by Voivodship Inspectorate of Environmental Protection (VIEP) in Gdańsk. We used data from stations located within the model effective domain (Figure 1). VIEP’s scope, methods of monitoring, and the manner of assessing the status of the Baltic Sea has been defined by the provisions of the Water Law Act and the monitoring program of marine waters adopted by the Council of Ministers, which implements the requirements of Article 11 of Directive of the European Parliament and of the Council 2008/56/EC of 17 June 2008, establishing a framework for community action in the field of marine environmental policy. Monitoring services performed by VIEP constitute the fulfilment of Poland’s obligations arising from the convention on the protection of the marine environment of the Baltic Sea area and, since 15 July 2014, of Article 11 of Directive 2008/56/EC of the European Parliament and establishing a framework for community action in the field of marine environmental policy.

2.5.2. NEMO-SCOBI

The source of the numerical results used in the manuscript is the data downloaded from the Marine Copernicus database (http://marine.copernicus.eu). This Baltic Sea biogeochemical reanalysis product provides a 25-year biogeochemical reanalysis for the Baltic Sea (1993–2017) using the ice-ocean model NEMO-Nordic (based on NEMO-3.6, Nucleus for European Modelling of the Ocean) coupled with the biogeochemical model SCOBI (Swedish Coastal and Ocean Biogeochemical model) together with LSEIK (local singular evolutive interpolated Kalman) data assimilation. All variables are available as daily means as well as monthly means and include nitrate, phosphate, ammonium, dissolved oxygen, and chlorophyll-a. The observation types used in the data assimilation are nitrate, phosphate, dissolved oxygen, and chlorophyll-a. The NEMO-SCOBI domain covers the whole Baltic Sea extended to the North Sea area.

2.6. Nutrient Spread Module

In order to model the spread of nutrients, the biochemical model was combined with a hydrological model describing the processes taking place in the Puck Commune’s watercourses and rivers in the model domain. For this purpose, we use (Equation (1)) in the Cartesian coordinate system without the source and sink term ($F_S$ = 0), where S is the concentration of nutrients and the rate of descent is $w_s=0$. The nutrient spread module describes the diffusion and advection of nutrients under the influence of sea currents. Data on the flow and concentration of nutrients are provided once a day at the land–water border from the SWAT model. The module consists of two partial second-order differential equations for concentrations of nutrients, e.g., phosphates and sum of nitrates and nitrites.

2.7. Pesticide Distribution

Models of pesticide distribution in surface waters have been reviewed with particular attention to the marine environment. As a result of the literature review, SWAT and FANTOM [13] models were selected as the basis for the development of the pesticide distribution module in the EcoPuckBay model. We developed a computational algorithm including the basic processes taking place in the sea. Sedimentation and degradation processes of pesticides have been introduced through a function of sources and sink $F_s=Q_p$ in the passive substance transport of Equation (1):

$$Q_p=w_s\frac{K_d{C}_s}{1+K_d{C}_s}\frac{\partial C_p}{\partial z}-k_{p,aq}{C}_p,$$

(2)

where $C_p$ is the concentration of pesticide and where empirical model parameters such as degradation factor $k_{p,aq}$ and adsorption factor $K_d$ have been calibrated in the SWAT hydrological model for surface water in watercourses. The concentration of the suspended matter $C_s$ plays a very important role in the adsorption of pesticides and was estimated on the basis of dead and live organic matter described in the biochemical model by appropriate state variables (describing plankton, detritus, etc.). Due to the large number of compounds considered as pesticides, the list of pesticides that are modelled in the Puck Bay area was narrowed down. The list includes pesticides currently used in agriculture. For selected compounds, surface water inflows of selected pesticides to the computational domain simulated with the use of SWAT hydrological model were estimated.

Pesticides are also taken from water by living organisms and accumulated in their bodies. However, since bioaccumulation has a small effect on the concentration of pollutants in water; it can be omitted in the first approximation when modelling pesticide transport in seawater.

2.8. Operational Mode

The hydrodynamic and biochemical parts of the EcoPuckBay model work in the operational version and are available on the website. Every 6 h, a 48-h forecast is generated based on the atmospheric forcing data provided by the UM-ICM weather model and data for the model domain’s borders. In addition, the model produces a forecast when new data is released from sources providing data for assimilation, e.g., SWAT, Modflow, 3D CEMBS, and satellite data (see Figure 2).

To connect with EcoPuckBay model data sets from various sources (characterised by very different data formats, spatial grid, and time coverage), scripts were created mainly in Python. They ensure autonomous and automatic operation of the WaterPUCK website with the ability to monitor the correctness of its operation by the administrator.

It is possible to create forecast maps of the Puck Bay for the following hydrodynamic parameters: temperature, salinity, currents, sea surface height, as well as biogeochemical parameters—chlorophyll-a concentration, nutrient concentrations (NO${}_3$, NH${}_4$, PO${}_4$, and SiO${}_3$), dissolved oxygen concentration, and active ingredients of pesticides. It is also possible to create graphs of temporal changes of these parameters at a given point as well as vertical cross sections. The results of simulations for a 48-h forecast are presented for the model with a vertical resolution of 115 m. The possibility of using the WaterPUCK service in operational mode can be helpful for authorities and policymakers to have access to scientific knowledge based on real environmental data as quickly as possible. In addition, prognostic information can be very useful for both shipping and broadly understood tourism sectors.

3. Results

This section gives some of the results obtained from the EcoPuckBay model and presents the Puck Bay’s coastal ecosystem functioning. The results presented below demonstrate that the model is operating correctly.

3.1. Statistical Comparison of EcoPuckBay Model With VIEP Monitoring Data and NEMO-SCOBI Model Results

To evaluate the quality of the results obtained from the EcoPuckBay model, a set of basic statistical measures were calculated, such as mean, Pearson correlation coefficient (r), root-mean-square-error (RMSE), and standard deviation (STD). The statistical comparison is presented below in the form of tables and modified Taylor diagrams [14]. Dissolved oxygen, nitrates, phosphates, and chlorophyll-a were selected for comparison.

3.1.1. VIEP Monitoring Data

By comparing the results of the EcoPuckBay model with the VIEP monitoring data, we see that the modelled variables concentrations are of the same order of magnitude as the environmental data. The highest uncertainty (about 30 per cent) is for nitrates; in other cases, the difference in averages does not exceed 10 per cent. In the case of temporal variability of the modelled variables, dissolved oxygen is best determined for which Pearson’s correlation coefficient is the highest ($r=0.71$). In

Table 3. Statistical comparison of EcoPuckBay model with Voivodship Inspectorate of Environmental Protection (VIEP) monitoring data.

	VIEP Mean	VIEP STD	EPB Mean	EPB STD	RMSE	r
O${}_2$ [mmol m${}^{-3}$]	315.30	43.06	337.50	35.97	31.03	0.71
NO${}_3$ [mmol m${}^{-3}$]	1.74	2.31	2.48	2.87	2.19	0.66
PO${}_4$ [mmol m${}^{-3}$]	0.35	0.23	0.32	0.21	0.21	0.56
CHL [mg m${}^{-3}$]	4.44	3.14	4.40	3.59	2.89	0.64

, we present the results of statistical analysis of the EcoPuckBay model in the context of VIEP environmental data.

Figure 3 presents in a graphical form information from Table 3 as a Taylor diagram (modified to show variables with different ranges of variation in one picture). As a rule, in the Taylor diagram, the closer the model point is to the reference point ($r=1,\mathrm{STD}=1,\mathrm{and}\mathrm{RMSE}=0$), the more accurately the model reproduces measurement data. From the diagram, it can be inferred that the model results of nitrates and chlorophyll-a are characterised by slightly greater inertia (greater range of variation) compared to environmental data while nitrates and dissolved oxygen vary less in the model than the measurement data. Root-mean-square-error (modified by dividing by the reference STD to compare all variables at one time) is less than 1.0 for all modelled variables.

3.1.2. NEMO-SCOBI Model Results

It should be noted that the EcoPuckBay and NEMO-SCOBI models have very different grid sizes. For this reason, for calculations, a sufficiently large number of EcoPuckBay cells were taken to cover the volume of a single NEMO-SCOBI cell. Analysing the statistical comparison of the EcoPuckBay model with the NEMO-SCOBI model, a high correlation for dissolved oxygen ($r=0.75$) can be observed, while the other variables have a Pearson correlation coefficient of about 0.6 (

Table 4. Statistical comparison of EcoPuckBay model with NEMO-SCOBI (Nucleus for European Modelling of the Ocean coupled with Swedish Coastal and Ocean Biogeochemical model) model data.

	NEMO Mean	NEMO STD	EPB Mean	EPB STD	RMSE	r
O${}_2$ [mmol m${}^{-3}$]	360.63	43.74	358.90	35.95	29.02	0.75
NO${}_3$ [mmol m${}^{-3}$]	4.68	5.18	5.70	3.82	4.20	0.60
PO${}_4$ [mmol m${}^{-3}$]	0.68	0.36	0.21	0.15	0.29	0.59
CHL [mg m${}^{-3}$]	6.37	5.87	3.51	2.40	4.74	0.63

).

A graphical representation of the comparison of both models is shown in Figure 4. The Taylor diagram was made in such a way that the reference point ($r=1,\mathrm{STD}=1, \mathrm{and}\mathrm{RMSE}=0$) corresponds to the NEMO-SCOBI model data. The RMSE is smaller than the reference data STD for all modelled variables (normalised RMSE in the diagram is smaller than 1.0). All presented variables have a smaller range of variations in the EcoPuckBay model compared to the NEMO-SCOBI model.

3.2. Land–Water Linkage

We coupled the EcoPuckBay model from the shore side with both the surface water runoff model based on SWAT and the groundwater flow model based on Modflow.

3.2.1. River Runoff

Figure 5 and Figure 6 present the average daily concentrations of nitrates and mineral phosphorus in 2016 for the waters of Reda, Gizdepka, Potok Błądzikowski, and Płutnica. The concentrations of the presented nutrients were determined using the coupling module with data provided by the SWAT hydrological model during 2016 and, as it can be observed in the following plots, varied greatly depending on the river.

The concentration of mineral phosphorus was highest in the waters of Reda and Płutnica, reaching up to 0.20–0.25 mg dm${}^{-3}$, while Gizdepka and Potok Błądzikowski carried waters with a mineral phosphorus concentration of about twice as low as 0.1–0.15 mg dm${}^{-3}$ (Figure 5). In addition, in the summer, it can be observed that higher concentrations are more frequent than in the other periods of 2016 and that this period has a much higher daily variability (peaks).

The time variability of the average modelled nitrates concentration is presented in Figure 6 and is qualitatively different for phosphates. In the first half of the year, concentrations remain fairly constant (with several clear dips for the Reda River). However, in the second part of the year, more dynamic changes dominate. The range of variability of mean concentrations for the modelled watercourses is in the range of 0–25 mg dm${}^{-3}$, the highest mean concentrations were found for the Reda River, and the lowest were found for the Płutnica River.

3.2.2. Groundwater Discharge

Modflow computations indicate that the rate of discharge of groundwater to Puck Bay (averaged monthly) varies from 1000 to ca. 1600 m${}^3$ h${}^{-1}$ (i.e., $7.3\times {10}^5$ to $11.7\times {10}^5$ m${}^3$ month${}^{-1}$). These results are consistent with previous studies on groundwater flow in the area [15]. The associated load of nitrate varies in the range of 452 to 1276 kg month${}^{-1}$. The results quoted here refer to the current conditions of land use and agricultural practices as obtained from the survey of Puck area. Modelling studies showed that changes in land use and agricultural practices have considerable influence on groundwater flow rates and nitrate loads at the land-sea interface [6].

3.3. Prognostic Model Parameters

We present the results for 2016 as a representative year for the last 5 years (differences in individual years do not exceed 10 per cent).

3.3.1. Dissolved Oxygen

Seasonal changes in water oxygenation are influenced by both climatic factors and the dynamics of primary production. The nature of these changes is typically seasonal (Figure 7). The maximum oxygen solubility occurs in winter when the water temperature is at its minimum. The highest monthly average of dissolved oxygen concentration occurred in February and was equal to 412.64 mmol dm${}^{-3}$.

Then, as the temperature increases, the solubility decreases, but in places where there is intensive primary production, increased levels of concentration are noticeable locally. The minimum concentration of dissolved oxygen in water occurs during the summer (specifically in August with a mean value of 313.37 mmol m${}^{-3}$) because of the lowest solubility due to the temperature maximum. The average annual dissolved oxygen concentration for the study area is about 358.06 mmol m${}^{-3}$, with a standard deviation of 37.96 mmol m${}^{-3}$.

3.3.2. Chlorophyll-a

The main source of energy for life on Earth is the Sun, the water environment is no exception. The scattering and reflective properties of organisms and molecules contained in the water column affect the depth of the euphotic layer, for which high values indicate the possibility of a primary production process in deeper parts of the sea. The rate of phytoplankton primary production in the growing season is usually very high. Owing to the short life span of these microscopic plants and the high productivity of the euphotic layer, phytoplankton is the main source of energy for other ecosystem elements. Some phytoplankton is consumed directly by herbivorous zooplankton, but a large part of the phytoplankton sinks to the bottom. Three phytoplankton abundance peaks are characteristic of the entire Puck Bay; they are determined by the phenologies of the various groups of algae, which are associated with the environmental requirements of these plants. Cyanobacterial blooms appear mostly in July, August, and September. The highest concentrations of chlorophyll-a are observed relatively close to the shore, where access to nutrients is greatest. Initially, no elevated chlorophyll-a concentration is observed in Puck Lagoon until June due to the geomorphological separation from the rest of the basin (Figure 8).

The annual average concentration of chlorophyll-a was 3.76 mg m${}^{-3}$, with a standard deviation of 2.82 mg m${}^{-3}$. The lowest monthly mean concentration of chlorophyll-a in the surface layer of the study area was determined in January and was equal to 0.44 mg m${}^{-3}$, while the highest monthly mean concentration of chlorophyll-a occurred in August and was equal to 8.39 mg m${}^{-3}$.

3.3.3. Nutrients

Nutrients such as nitrogen and phosphorus, besides atmospheric deposition, can enter the marine environment from two types of sources: point and nonpoint sources. Point sources are usually municipal wastewater treatment plants, while nonpoint sources include surface runoff from agricultural fields. Therefore, nitrogen and phosphorus are the main limiting factors of biological production and, therefore, it is important to know both the ratio of these two elements and its variability over time. The nutrient concentration in the Puck Bay shows seasonal variability. The maximum concentrations are usually recorded in winter and early spring, before the start of the vegetation period, which starts immediately after the ice has descended. Knowledge of concentration and nutrient load is essential for assessing water quality and for understanding ecosystem functioning. The total nutrient load entering the Bay of Puck is more important than their concentration in watercourses due to advection and diffusion. The load discharged by the examined streams depends on their outflow, which varies significantly throughout the year.

Nitrates

The highest concentrations of nitrate–nitrogen were found in winter and early spring before the beginning of the growing season, and the lowest was found in the middle of the growing season (Figure 9).

The average annual nitrate concentration in the surface layer was 5.06 mmol m${}^{-3}$, with a standard deviation of 3.90 mmol m${}^{-3}$. The highest average monthly nitrate concentration was determined in January (11.04 mmol m${}^{-3}$), while the lowest nitrate concentration occurred in August (1.29 mmol m${}^{-3}$).

Phosphates

Phosphate concentrations in the Puck Bay waters also show seasonal variability, related to the size of loads supplied by surface waters and the depletion of these compounds during the vegetation period as a result of phytoplankton consumption. A clear decrease in the concentration of nitrogen compounds was observed in the summer months, and in early autumn, it increased again (Figure 10).

The annual mean concentration of phosphates in the surface layer was 0.21 mmol m${}^{-3}$, with a standard deviation of 0.14 mmol m${}^{-3}$. The highest monthly average phosphate concentration occurred in December (0.51 mmol m${}^{-3}$), while the lowest phosphate concentration occurred in August (0.04 mmol m${}^{-3}$).

3.4. Nutrient Spread

Figure 11 shows the intensity and character of the spread of nutrients and pollutants supplied to the Baltic Sea waters by rivers flowing within the model domain. Due to the relatively even but dense distribution of rivers along the west coast of the domain (from Gdańsk-Vistula to Puck/Władysławowo, where Płutnica flows), there is a high probability that, after just 1–2 days, the pollution entering the rivers as a result of mixing and currents will be spilt along the whole west coast. On day 3, nutrients can be transported to the coast of Hel Peninsula and will spill along its coastline in the following days. At the same time, with the rapid transport of the pollutants along (and in close proximity) to the shore, they spill over into the central areas of the Puck Lagoon, although it usually takes at least 3–4 days for the pollutants to reach the areas of the Puck Bay lying at an equal distance between the west coast and the Hel Peninsula, although this is not very likely to be the case. After 5 days, there is a nonzero probability that a nutrient or other pollutant that has entered the sea from rivers in the domain will reach any part of the Puck Bay. From day 6, it is visible that, in the central part of the bay, there is an increased probability of nutrients appearing up to 50–60 per cent in the central part of the bay and over 90 per cent in all waters in close proximity to the coastline. There is also a strong expansion of nutrients coming from the Vistula as a result of the amount of water pumped by the Vistula and the sometimes favourable physical conditions (currents causing water to be pumped inside the Puck Bay).

Analysing the distribution of currents in individual months of the year 2016 (see Figure 12), we can observe two possible scenarios: the first one when the water is pushed into the bay (i.e., October) and the second when the water is drained from the bay (i.e., August and December). In the first case, sea surface height on the western side of the domain is higher than on the eastern side (see [1]). In the second case, the sea surface height is relatively lower on the western side of the domain than on the eastern side.

The north direction of the currents along the west border of the domain is noticeable. The consequences of this state are visible on nutrient distribution maps. Namely, biogenic substances flowing with river waters into the Baltic Sea are often transported along the coast to the north and, only after some time, they spread into the central waters of the Puck Bay. Pollutants, therefore, tend to spread rapidly in the close coastal zone and to expand more slowly into deeper waters far from the coast. One can also see a great influence of biogenic substances coming from the Vistula, which together with the currents according to the scenario when water is pumped into the bay are transferred in the northern and northwestern direction, thus feeding the Puck Lagoon with pollutants and biogenic substances not coming from the rivers directly entering the Puck Lagoon.

Figure 13 presents how many days must elapse from the beginning of the distribution to the moment when a nutrient with a given probability or more will appear at a specific point in the domain. For example, “$p>10\%$” shows a situation in which, with a probability of 10 per cent or more (i.e., at least once in ten cases), the nutrient from rivers reaches a given point on the bay. Results show that 9 days is enough to reach every place (even the inner part) in the bay with a probability higher than 10 per cent. The increased probability thresholds are intended to show the regions along the west coast of the domain where polluted waters are most likely to reach after only a few days due to northward currents. Spread of pollutants into the central parts of the bay occurs either after a long time or less frequently—i.e., as can be seen in the case “$p>50\%$”, there are places (marked in white) where the nutrient from the rivers did not appear even after 9 days in half the cases.

It should be emphasised that biological processes are not taken into account in the nutrient spread module (expression $F_S=0$ in Equation (1)). Such an approach results in the maximum possible contamination, as the compound could be found to be decomposed and undetectable under real conditions (e.g., in the middle of the bay).

3.5. Pesticides

The concentrations of detected pesticides in soil samples (chlorpyrifos-ethyl, antrachinon, glyphosate, and AMPA (aminomethylphosphonic acid)), measured within the other stage of the WaterPUCK project, ranged from 0.05 to 0.35 mg kg${}^{-1}$. The highest concentrations were obtained in August in water from drainage ditches for metazachlor (2.0 mg dm${}^{-3}$), glyphosate (4.7 mg dm${}^{-3}$), and AMPA (2.0 mg dm${}^{-3}$). In the surface water samples and sediment samples taken from the Puck Bay, none of the 309 substances from the pesticides group were found [16]. However, 5 pesticides of agricultural origin, which are glyphosate, diflufenican, metazachlor, chlorpyrifos, and antrachinon, are modelled in the EcoPuckBay model. In the online portal, it is possible to display spatiotemporal distributions for the concentration of glyphosate, which is the most commonly used pesticide in the studied area. The distributions of other modelled pesticides are very similar, which is why we present only one example on the webpage.

Figure 14 shows a screenshot of the web portal showing the spatial distribution map for the active substance of the pesticide being modelled. It can be seen that the concentrations are so low that it is not possible to detect this compound in the waters of the Gulf of Puck by easily accessible methods.

3.6. Web Portal

Nowadays, quick access to expert knowledge is highly valued, especially in the context of decision-making by the authorities. In order to meet these requirements, we developed this internet service. This is one of the key services created within the framework of the project on which the results of all models included in the WaterPUCK Integrated Information and Prediction Service are available. This website operates dynamically in the operational mode, allowing for visualisation of forecast maps, time series, and vertical sections.

In this paper, we present the website which allows to generate maps of selected biochemical parameters from the EcoPuckBay model in the area of the Puck Bay and the western part of the Gulf of Gdańsk. It is now possible to generate raster maps (Figure 15a) for individual depths, which represent the next vertical level of the model. In addition, it is possible to create a time series (Figure 15b) for the set periods in the selected location (after prior determination or indication of the desired latitude and longitude) as well as W-E and/or S-N cross sections, allowing for analysis of the variability of the parameter state in the entire water column (Figure 15c).

An important advantage of the portal is relatively fast access to archival model data and 48-h forecasting data. Using the calendar, any date from the modelled period is selected in an intuitive way. This service can be accessed from the project website (www.waterpuck.pl) after choosing “Products” from the navigation menu and by selecting “EcoPuckBay Biochemical Model of the Puck Bay”.

4. Discussion

4.1. Prognostic Model Parameters

The Puck Bay, and especially its inner part—the Puck Lagoon—is a very dynamic area with relatively low inertia. Evaluation of the hydrodynamic part of the model [1] showed that, in this part, the annual differences in salinity and temperature are the highest. There are several models describing processes occurring in the marine environment in use such as HIROMB (High-Resolution Operational Model for the Baltic Sea), RCO (Rossby Centre Ocean Model), NEMO (Nucleus for European Modeling of the Ocean) and biogeochemical model SCOBI (Swedish Coastal and Ocean BIogeochemical model), HIRLAM (High-Resolution Limited Area Model), or the ecosystem model of the Baltic Sea BALTSEM. At the end of 2015, the SatBałtyk System has been initiated in the IO PAN (Institute of Oceanology of the Polish Academy of Sciences). It should be noted that, for this area, this is the first model with such high resolution that also is combined with terrestrial water flow models (such as the SWAT and Modflow). Also, earlier studies conducted by [17] show that the studied area is extremely difficult to model due to the complexity and intensity of processes occurring in its ecosystem.

The annual averages are of the same order of magnitude for all data sets (

Table 5. Annual means of the selected variables.

	O${}_2$ [mmol m${}^{-3}$]	NO${}_3$ [mmol m${}^{-3}$]	PO${}_4$ [mmol m${}^{-3}$]	CHL [mg m${}^{-3}$]
EPB mean	358.06	5.06	0.21	3.76
NEMO mean	360.63	4.68	0.68	6.37
VIEP mean	315.30	1.74	0.35	4.44

). For dissolved oxygen, the average of the EcoPuckBay model is between the value determined from the NEMO-SCOBI model and from the VIEP measurement data. For nitrates, the average of the EcoPuckBay model is higher than that of the NEMO-SCOBI and the VIEP measurements. For phosphates and chlorophyll-a, on the other hand, the EcoPuckBay model has a lower average compared to the averages determined from the NEMO-SCOBI model and VIEP environmental data. It should be noted that the data from the NEMO-SCOBI model as well as the VIEP measurement data have a much lower spatial and temporal coverage than the data from the EcoPuckBay model.

A very characteristic picture for the spatial and temporal changes for all the examined variables in the surface layer is the very strong separation of the outer part of Puck Bay from the Puck Lagoon. This division is most visible when nutrients from the Vistula River runoff move with the strong eastern winds deeper into the Puck Bay. In such a situation, there is a significant (visible on maps) difference between the concentrations of nutrients in Puck Lagoon and those in the outer Puck Bay in a short time.

The vertical distributions of modelled variables are mainly determined by hydrodynamic processes occurring in the Puck Bay. As a rule, heated and less salty water on the surface is characterised by a high concentration of chlorophyll-a and a lower concentration of nutrients (which were used for primary production). Furthermore, during the growing season, this condition intensifies and the surface waters are clearly separated from the bottom waters rich in nutrients, where decomposition and mineralisation occur.

From the results obtained, it can be concluded that the overall level of impact of farms from the Puck Commune on the environment, including the quality of the waters of the Puck Bay, is relatively low. It should be noted, however, that among the studied farms, there were some for which the environmental pressure due to surveys can possibly have an unusually huge impact.

4.2. Nutrient Spread and Pesticides Distribution

The methodology used in the “EcoPuckBay—Nutrient spread module” assumes visualisation of the distribution of biogenic substances. In real environmental conditions, the presented results of the probability distribution may be different due to the fact that the biogenic substances carried by sea waters may be used by phytoplankton for growth in the primary production process and that, at the same time, their concentrations may increase as a result of the processes of water masses and chemical processes, i.e., remineralization in the water column resuspension from bottom sediments and secretion by zooplankton. The model shows the potential maximum range to which nutrients can spread and their concentrations if no consumption or growth were to occur as a result of appropriate environmental conditions.

As a result of the conducted studies, we come to the conclusion that active substances contained in pesticides used in the agricultural sector in the Puck Commune are present in the marine environment of the Puck Bay in concentrations so small that they are undetectable. Moreover, we come to the same conclusions by conducting numerical research.

4.3. Web Portal

The WaterPUCK’s web portal is an excellent tool that allows to intuitively track changes in both archival and forecasting of selected biochemical variables. In particular, it allows the presentation of data in specific depth layers and vertical sections (see Figure 16, Figure 17 and Figure 18).

Using the data presentation tools provided as part of the WaterPUCK service, it is possible to check the status of the waters of the Puck Bay. For example, in Figure 16, we see how the concentration of chlorophyll-a decreases with depth, while in Figure 17, we see a higher concentration of nitrates at the bottom. Such a picture reflects the structure of the primary production process taking place, during which nitrates are consumed in the euphotic layer.

Furthermore, it is possible to study the spread and accumulation of pesticides in various parts of the Puck Bay (Figure 18).

5. Conclusions

High-resolution mapping of the state of the marine environment and its ecosystem was, is, and will probably be the object of interest not only for researchers but also for policymakers and authorities. The development of the WaterPUCK service meets these requirements and offers a product that makes it possible to monitor and forecast the state of the marine environment of the Puck Bay in a resolution that has not been offered by any other models so far. The key features of the WaterPUCK system are the use of data from land water models (surface water (SWAT) and groundwater (Modflow)) fed with data from farmers’ surveys and environmental measurements. This approach to modelling, especially in the context of modelling coastal waters, is very desirable because, in such a small and dynamically changing water body, correct mapping is only possible taking into account the vast majority of inputs and processes occurring.

The EcoPuckBay model is a useful tool allowing to control the state of the environment of the Puck Bay in terms of hydrodynamic and biochemical parameters, including pollution with nutrients and pesticides based on numerical studies and indirectly through the assimilation of satellite data. It allows tracking changes occurring in the studied area to improve the condition and protection of the ecosystem as well as predicting its response to threats arising from natural causes, independent of human activity as well as threats caused by human activities.

This article presents the validation of the three-dimensional ecohydrodynamic model of the Puck Bay named EcoPuckBay (specifically its ecosystem part). It is possible, via the website www.waterpuck.pl, to access archival and forecast data. The biochemical part presents the results for seven variables (nitrate, phosphate, ammonia, silicate, chlorophyll-a, dissolved oxygen, and active ingredient of pesticide). Presentation of spatial distributions for any depth, time variability for a selected point, and a vertical cross section between any two points within the domain set by the user is possible. Comparison of the results of the EcoPuckBay model with real measurement data as well as with the NEMO-SCOBI model data indicate the correct parametrization of the EcoPuckBay model.

The model contains a nutrient spread module with which it is possible to track the fate of agricultural nutrients in the coastal environment of the Puck Bay. The presented approach to the description of pollution distribution may be of particular importance in anticipation of events resulting in the occurrence of high-risk states or even ecological disasters by selecting potentially dangerous areas or by selecting safe areas.

The EcoPuckBay model can certainly be further improved. The element that is probably most likely to improve the quality of model performance is to introduce the Vistula River into the model based on actual hydrological conditions, as the studies have shown that its impact on the water environment of the Puck Bay is the greatest.