Integrated Hydrodynamic and Atmospheric Modeling of Polyfytos Lake for Substance Dispersion Using Delft3D and Weather Research and Forecasting Model

Vassilios Papaioannou^*, Damianos Florin Mantsis, Christos Ge Anagnostopoulos, Konstantinos Vlachos, Anastasia Moumtzidou, Ilias Gialampoukidis, Stefanos Vrochidis, Ioannis Kompatsiaris

Information Technologies Institute, Centre for Research and Technology Hellas, Thessaloniki, Greece.

DOI: 10.4236/ojce.2025.152013                

Abstract

This study investigates the hydrodynamic and environmental behavior of Polyfytos Lake in Greece using the three-dimensional modeling tools developed in Delft for water and environmental systems (Delft3D), including the particle tracking module (Delft3D-PART) for simulating the transport and dispersion of pollutants. The primary objective was to simulate water movement and assess the dispersion of a hypothetical pollutant, such as oil, under realistic atmospheric conditions. Meteorological forcing was provided by high-resolution data derived from the Weather Research and Forecasting (WRF) model, ensuring an accurate representation of wind, temperature, and precipitation patterns across the lake surface. The integration of WRF data into Delft3D enabled dynamic and time-dependent boundary conditions that significantly enhanced the reliability of the simulation. A hypothetical heavy oil spill was introduced, with specific properties designed to emulate the behavior and transport of an oil-like substance. The Delft3D-PART module was employed to model the particle-based transport of the pollutant, tracking its movement, dispersion, and interactions with the lake’s hydrodynamic conditions over time. Its distribution was monitored to evaluate potential environmental impacts and inform response strategies. The study demonstrates the utility of Delft3D as a decision-support tool for environmental risk assessment and highlights the importance of incorporating realistic meteorological forcing into hydrodynamic models for effective management of complex inland water bodies.

Keywords

Delft3D-FLOW, Delft3D-PART, WRF Model, Oil Spill, Polyfytos Lake, Simulation

1. Introduction

Lakes are vital ecosystems that support water supply, biodiversity, and diverse human activities, including agriculture, recreation, and energy production [1] [2]. Understanding water circulation and pollutant dispersion within these systems is essential for environmental protection and informed decision-making. To achieve this, researchers often rely on computer models that simulate hydrodynamic processes and environmental changes over time [3].

One widely used tool is Delft3D, a modeling software suite developed by Deltares [4], which simulates water flow, sediment transport, and pollutant dispersion in aquatic environments such as rivers, lakes, and coastal zones [5]. Its applications span various settings, including subtropical estuaries like the Vitória Island Estuarine System (VIES) in Brazil, where computational fluid dynamics coupled with Reynolds-averaged Navier-Stokes equations have effectively modeled water quality stresses under significant anthropogenic pressure [6]. Hydrodynamic models are also increasingly applied to urban lakes, which provide critical ecosystem services such as flood control, nature conservation, and recreation. For example, Lake Créteil in France—a small, shallow urban lake—was studied using the Delft3D-FLOW module, calibrated with high-resolution, high-frequency measurements to simulate temperature dynamics, nutrient diffusion, and pollutant resuspension [7].

Accurate Delft3D simulations require detailed meteorological input, known as atmospheric forcing, which includes wind speed and direction, air temperature, and precipitation. The Weather Research and Forecasting (WRF) model is commonly used to generate such high-resolution, region-specific weather data [8]. However, because WRF and Delft3D are not directly compatible, WRF outputs must be preprocessed and converted into Delft3D-compatible formats. While coupling WRF and Delft3D has been explored in prior studies for coastal [9] or estuarine [10] applications, its implementation in an inland, thermally stratified reservoir such as Polyfytos Lake remains relatively underexplored. Moreover, the development of a reusable Python-based preprocessing pipeline enhances replicability and operational utility, particularly in regions lacking integrated meteorological-hydrodynamic systems.

Recent urban meteorology research highlights the growing importance of high-resolution data for sustainable decision-making. For instance, spatio-temporal maps of near-surface urban air temperature (T_a) are crucial for assessing urban heat exposure and climate impacts. A study in Warsaw, Poland, combined crowd-sourced Netatmo weather station data with machine learning and traditional interpolation to improve daily T_a predictions. The WRF model also demonstrated that machine learning offers an efficient alternative for identifying urban heat-vulnerable areas amid climate change [11]. Beyond urban climate, WRF-Chem coupled with the MEGANv2.1 emission module was used to estimate biogenic volatile organic compound emissions in the Colombian Andes, illustrating WRF’s versatility in ecological modeling [12]. Furthermore, spatial resolution plays a key role in regional strategic environmental assessments (R-SEA), affecting atmospheric prediction outcomes when varying emission data scales [13]. In watershed management, comparisons of observed and WRF-simulated meteorological datasets showed that weather data variations can significantly impact nitrogen fate and transport modeling using tools like SWAT [14].

This study focuses on Polyfytos Lake, Greece’s largest artificial reservoir, located in Western Macedonia. Created in 1975 by damming the Aliakmon River, the lake covers approximately 74 km². It is vital for hydroelectric power generation, irrigation, and—since 2003—supplying drinking water to Thessaloniki, Greece’s second-largest city with around 1,050,000 residents. Approximately 145,000 m³ of surface water is extracted daily to meet the city’s demand [15].

The objectives of this study are to:

1) Simulate water circulation in Polyfytos Lake and examine contaminant dispersion (e.g., an oil spill) under varying conditions using the Delft3D-PART module.

2) Integrate realistic weather patterns by incorporating WRF-generated atmospheric data.

3) Develop and demonstrate a Python-based workflow for converting WRF outputs into Delft3D-compatible input files.

This replicable workflow enables researchers and practitioners to apply the approach to other lakes or geographic areas, thereby improving the realism and accuracy of model simulations to support enhanced environmental planning, risk assessment, and emergency response strategies.

2. Model Specification and Numerical Procedure

This section provides a detailed overview of the coupled hydrodynamic and atmospheric modeling approach applied to Polyfytos Lake. We begin by introducing the study area, highlighting its key geographical and environmental characteristics. Next, we present the governing equations that underpin the hydrodynamic simulations specific to this region. Following this, the development and validation of the hydrodynamic computational grid are described to ensure the accuracy and reliability of the simulation results. Finally, we outline the setup and validation of the Weather Research and Forecasting (WRF) model, whose outputs serve as the atmospheric forcing data integrated into the hydrodynamic model.

2.1. Case Study Area

Polyfytos Lake, located in the Kozani regional unit of northern Greece, is the country’s largest artificial reservoir. It was created in 1974 with the construction of the Polyfytos Dam on the Aliakmon River, primarily to support hydroelectric power generation [16]. The lake covers a surface area of approximately 74 km² and has a maximum storage capacity of 1100 million cubic meters. The associated hydroelectric power station has a design capacity of 375 MW. The surrounding region is noted for its distinctive geographical features and ecological richness (Figure 1).

Figure 1. Map of the Polyfytos Lake study area, including the location of nearby meteorological stations used for atmospheric data input.

Around 145,000 m³ of surface water is withdrawn daily to supply the city’s Drinking Water Treatment Plant [17]. In addition to its utilitarian functions, the lake supports diverse socio-economic and environmental activities. The Environmental Education Center of Lake Polyfytos highlights the importance of conserving the region’s natural habitats and biodiversity. The area is also home to significant infrastructure, including the Polyfytos Bridge, which spans the lake and plays a vital role in regional connectivity and transportation. The bridge and its surrounding environment have been the focus of multiple studies addressing structural integrity, environmental impacts, and ecosystem preservation [18].

2.2. Governing Equations

The governing three-dimensional equations describing free-surface flows can be derived from the Navier-Stokes equations after averaging over turbulence time-scales (Reynolds averaged Navier-Stokes equations). Such equations express the physical principle of conservation of volume, mass and momentum [4]. In this section we describe the shallow water equations, for which the depth is assumed to be much smaller than the horizontal length scales of flow and bathymetry. Under the further assumption that the vertical accelerations are small compared to the horizontal ones the shallow water assumption is valid. This means that the vertical momentum equation is reduced to the hydrostatic pressure relation.

$\frac{\partial p}{\partial z}=-{\rho }_o\cdot g$ (1)

where: p = pressure, ρ_ο = the reference water density and g = gravity acceleration.

To facilitate the representation of flows over variable bathymetry, the vertical coordinate is transformed into a terrain-following σ-coordinate system, defined such that σ = 0 at the free surface and σ = −1 at the bottom. The three-dimensional hydrostatic shallow water equations, which for convenience of presentation are given in Cartesian rectangular coordinates in the horizontal and σ-coordinates in the vertically.

The horizontal momentum equations are:

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\frac{w}{h}\frac{\partial u}{\partial \sigma }=-\frac{1}{{\rho }_o}\left(\frac{\partial p}{\partial x}+\frac{\partial \sigma }{\partial x}\frac{\partial p}{\partial \sigma }\right)+fv+F_u+\frac{1}{h^2}\frac{\partial }{\partial \sigma }\left(v_V\frac{\partial u}{\partial \sigma }\right)$ (2)

$\frac{\partial \upsilon }{\partial t}+u\frac{\partial \upsilon }{\partial x}+\upsilon \frac{\partial \upsilon }{\partial y}+\frac{w}{h}\frac{\partial \upsilon }{\partial \sigma }=-\frac{1}{{\rho }_o}\left(\frac{\partial p}{\partial y}+\frac{\partial \sigma }{\partial y}\frac{\partial p}{\partial \sigma }\right)-fu+F_{\upsilon }+\frac{1}{h^2}\frac{\partial }{\partial \sigma }\left(v_V\frac{\partial v}{\partial \sigma }\right)$ (3)

The continuity equation (mass conservation) is:

$\frac{\partial \zeta }{\partial t}+\frac{\partial \left(hu\right)}{\partial x}+\frac{\partial \left(hv\right)}{\partial y}+\frac{\partial w}{\partial \sigma }=h\left(q_{in}-q_{out}\right)+P-E$ (4)

where: u, υ, w = the velocity components in the x, y, and σ directions respectively, h = the total water depth, f = the Coriolis parameter accounting for Earth’s rotation, F_u, F_υ = horizontal viscosity terms, υ_V = the vertical eddy viscosity, q_in and q_out = local water source and sink terms per unit volume (e.g., river inflows, withdrawals), P and E denote precipitation and evaporation fluxes at the surface.

The vertical velocities w in the σ-coordinate system are computed from the continuity equation:

$\frac{\partial \zeta }{\partial t}+\frac{\partial \left(hu\right)}{\partial x}+\frac{\partial \left(hv\right)}{\partial y}+\frac{\partial w}{\partial \sigma }=Q$ (5)

where: Q represents the net contributions per unit area from discharge, withdrawal, precipitation, and evaporation, defined as:

$Q=h{\displaystyle {\int }_{-1}^0\left(q_{in}-q_{out}\right)\text{d}\sigma }+P-E$ (6)

The horizontal viscosity terms F_u and F_υ account for turbulent momentum fluxes and are expressed as:

$F_u=\left(\frac{\partial }{\partial x}+\frac{\partial \sigma }{\partial x}\frac{\partial }{\partial \sigma }\right){\tau }_{xx}+\left(\frac{\partial }{\partial y}+\frac{\partial \sigma }{\partial y}\frac{\partial }{\partial \sigma }\right){\tau }_{xy}$ (7)

$F_{\upsilon }=\left(\frac{\partial }{\partial x}+\frac{\partial \sigma }{\partial x}\frac{\partial }{\partial \sigma }\right){\tau }_{xy}+\left(\frac{\partial }{\partial y}+\frac{\partial \sigma }{\partial y}\frac{\partial }{\partial \sigma }\right){\tau }_{yy}$ (8)

where: the Reynolds stresses τ_xx, τ_xy and τ_yy are parameterized using the Boussinesq approximation:

${\tau }_{xx}=2v_H\left(\frac{\partial u}{\partial x}+\frac{\partial \sigma }{\partial x}\frac{\partial u}{\partial \sigma }\right)$ (9)

${\tau }_{yy}=2v_H\left(\frac{\partial v}{\partial y}+\frac{\partial \sigma }{\partial y}\frac{\partial v}{\partial \sigma }\right)$ (10)

${\tau }_{xy}=v_H\left(\frac{\partial u}{\partial y}+\frac{\partial \sigma }{\partial y}\frac{\partial u}{\partial \sigma }+\frac{\partial v}{\partial x}+\frac{\partial \sigma }{\partial x}\frac{\partial v}{\partial \sigma }\right)$ (11)

where: v_H = horizontal eddy viscosity coefficient.

2.3. Grid Formation and Validation

Grid formation is a fundamental step in setting up hydrodynamic and environmental simulations in Delft3D for aquatic systems such as lakes, rivers, and coastal zones. The computational grid defines the spatial resolution and geometric structure over which key physical processes—such as water flow, dispersion, and pollutant transport—are calculated. For the Polyfytos Lake simulations, constructing a detailed and well-structured grid was essential to accurately capture the lake’s bathymetry, hydrodynamic behavior, and interaction with atmospheric forcing.

The first stage involved defining the spatial domain of the lake (Figure 2). For Polyfytos Lake, this required delineating the grid boundaries to encompass the full extent of the lake’s surface area. Delft3D’s Open Land Boundaries were employed to outline the lake’s natural shape and depth contours. These boundaries were supplemented by Open Samples to account for bathymetric variations within the domain.

Figure 2. RGFGRID interface in Delft3D.

The grid boundaries were defined using a Cartesian coordinate system and formatted for compatibility with Delft3D’s RGFGRID (version 6.03.00.66818), as illustrated in Figure 3(a). This grid file format includes metadata such as the file version and creation date, followed by a sequence of boundary points recorded in xyz coordinates. Specifically,

• x and y represent spatial positions (longitude and latitude, Figure 3(a)), and

• z denotes elevation or water depth at that location (Figure 3(b)).

• This approach ensures that the model accurately reflects the lake’s geometry, which is essential for simulating realistic hydrodynamic conditions. It is important to note that the Delft3D-PART module, which simulates particle transport, operates only in a Cartesian coordinate system.

Figure 3. Left: Land boundary definition for Polyfytos Lake, outlining the spatial limits of the simulation domain. Right: Bathymetry data for the lake, depicting depth variations used to model hydrodynamic conditions.

Once the computational grid was established, the next critical step was to define the boundary conditions required for simulation. These include parameters such as inflow and outflow discharges, atmospheric forcing variables (e.g., wind, temperature, precipitation), and water quality characteristics. Accurate and properly formatted boundary conditions are essential for producing realistic and reliable model outputs.

To streamline this process, the QUICKIN module within Delft3D was utilized. QUICKIN facilitates efficient setup and assignment of boundary conditions, ensuring input data is correctly formatted and aligned with the grid structure (Figure 4). This integration maintains spatial and temporal consistency across the model domain.

To assess the impact of grid resolution on temperature simulation accuracy in Polyfytos Lake, three spatial resolutions were tested: 200 m × 200 m, 100 m × 100 m, and 50 m × 50 m. Each resolution provided a different level of spatial detail, influencing the model’s ability to capture thermal gradients and local hydrodynamic phenomena.

Figure 4. Correctly formatted boundary condition data integrated with the Delft3D grid.

The model grid consists of 236 points in the M-direction and 233 points in the N-direction, striking a balance between computational efficiency and spatial detail. Vertically, the model is discretized into 14 layers with non-uniform thickness designed to enhance resolution near the surface where temperature gradients are more pronounced. Specifically:

• Layers 1 to 6: 0.2% thickness each,

• Layer 7: 1% thickness,

• Layer 8: 2% thickness,

• Layer 9: 5% thickness,

• Layer 10: 10% thickness,

• Layers 11 to 14: 20% thickness each.

This layered structure ensures high resolution in the thermally dynamic surface region while optimizing computational resources at greater depths.

Temperature profiles were analyzed at a representative point within the domain—coordinates M136, N123—corresponding to a critical observation site labeled OilSpill3out (Figure 5). This location was strategically selected to evaluate temperature distribution across the transition zone from deep to shallow waters, facilitating comparison across different grid resolutions (Figure 6).

Simulations covered the period from 1 July 2018, 00:00 UTC to 7 July 2018, 23:00 UTC, a timeframe characterized by prolonged heat and dry weather conditions. Temperature profiles were extracted for 6 July 2018, chosen to illustrate the effects of sustained elevated temperatures on thermal stratification.

Figure 5. Selected OilSpill3out monitoring point within Polyfytos Lake used for temperature profile analysis across different grid resolutions.

Figure 6. Temperature profiles at the OilSpill3out monitoring point for the three grid resolutions (200 m × 200 m, 100 m × 100 m, and 50 m × 50 m) in Polyfytos Lake.

The results demonstrate clear thermal stratification typical of medium-depth lakes:

• Surface temperatures started at approximately 24.9˚C, consistent with expected seasonal and atmospheric conditions [19].

• The thermocline, between 10 m and 30 m depth, exhibited a sharp temperature decrease from around 21.2˚C to 6˚C, indicating strong thermal gradients.

• Below 30 m, temperatures stabilized between 4˚C and 6˚C in the bottom layer, reflecting relatively constant cold conditions.

These temperature profiles provide vital insights into heat distribution and thermal dynamics, which are crucial for water quality, ecological health, and pollutant transport modeling. The profiles generated from the three tested grid resolutions differed by less than 1.0%, indicating robust reproduction of thermal stratification across all grids. Considering the trade-off between computational demand and accuracy, the 200 m × 200 m grid resolution was selected as optimal for the large-scale study, offering reliable results without compromising model precision.

2.4. WRF Setup and Validation

To improve the accuracy of Delft3D hydrodynamic simulations for Polyfytos Lake, atmospheric conditions are supplied by the Weather Research and Forecasting (WRF) model, version 4.3 [8]. The WRF simulations employ multiple nested domains with progressively finer resolutions to capture atmospheric processes at varying scales:

• d01: Large-scale domain with 25 km horizontal resolution,

• d02: Nested intermediate domain with 5 km resolution,

• d03: Inner domain with 1 km resolution,

• d04: Large Eddy Simulation (LES) domain with 200 m resolution.

The WRF model uses advanced physical parameterizations to represent sub-grid scale processes, including:

• The two-moment Morrison scheme for microphysics [20],

• The RRTMG scheme for longwave and shortwave radiation [21],

• The Betts-Miller-Janjic (BMJ) scheme for convection [22],

• The Mellor-Yamada-Janjic (MYJ) scheme for the Planetary Boundary Layer (PBL) [23],

• The Noah Land Model for land surface processes.

The 200 m resolution WRF domain (d04) was validated using observational data from two local meteorological stations—Kesaria and Velventos—situated near Polyfytos Lake. The meteorological data, provided by the National Observatory of Athens [24], include 10-minute averaged wind speed and dominant wind direction records spanning 1 July 2018 00UTC to 9 July 2018 00UTC. Table 1 summarizes the locations and elevations of these stations:

Table 1. Meteorological stations adjacent to Lake Polyfytos.

Station Code	Name	Longitude (˚)	Latitude (˚)
344	Velventos	22.056100	40.255576
449	Kesaria	21.859620	40.183720

The comparison covers 11,520 minutes (~8 days). Wind direction data, expressed as dominant directions, are shown every 960 minutes (~16 hours). The WRF-simulated wind speed and direction time series closely match observations from both stations, confirming the model’s capability to reliably reproduce local meteorological conditions (Figure 7(a) & Figure 7(b)).

Two custom Python scripts were developed to support the data processing workflow. The first script, wrf_degrees2cartesian.py (available on GitHub), converts the WRF model output from geographic (latitude–longitude) to Cartesian coordinates, enabling integration with the hydrodynamic model grid. The second script, extract_values_wrf.py, extracts wind speed data from WRF NetCDF output files. It calculates the mean wind speed at a specified location near the center of Polyfytos Lake and formats the data for direct comparison with observational records from the nearby meteorological station.

(a)

(b)

Figure 7. Wind speed and direction comparison between local stations (a) Kesaria station and (b) Velventos station and WRF model output.

While the validation of WRF meteorological inputs demonstrated strong agreement with observational data, further validation steps are planned to evaluate the hydrodynamic model against measured water levels and current velocities. Where available, comparisons with satellite or drone-based oil slick observations will also be integrated to benchmark the model’s ability to capture spatial spill patterns. These additions will enhance the model’s credibility for predictive and planning purposes.

3. Methodology

This section presents the comprehensive approach used to model the hydrodynamics and simulate oil spill behavior in Polyfytos Lake. The methodology integrates the Delft3D modeling system with meteorological forcing from the WRF model to realistically represent water flow, mixing processes, and the transport and fate of oil under varying environmental conditions. Detailed descriptions of the hydrodynamic setup, boundary conditions, and particle tracking techniques are provided to outline the simulation framework and parameterizations applied throughout the study.

3.1. Delft3D-FLOW Setup

A) Domain: The hydrodynamic modeling of Polyfytos Lake was conducted using the Delft3D-FLOW module, operating in a three-dimensional baroclinic mode. The model domain was discretized on a rectangular Cartesian grid, previously defined with 236 and 233 grid points in the M- and N-directions, respectively. The vertical structure consisted of 14 sigma layers, with thickness percentages described in Section 2.3.

B) Time frame: The simulation was conducted over the period from 01 July 2018 00:00 UTC to 08 July 2018 00:00 UTC (~7days), capturing a full week of meteorological and hydrodynamic variability relevant for the subsequent dispersion modeling. The computational time step was set at 1 min, selected based on a Courant–Friedrichs–Lewy (CFL) condition analysis to ensure numerical stability of the hydrodynamic solver. The CFL condition [4], defined as:

$CFL=\frac{U_{flow}^{max}\cdot \Delta t}{\Delta x}$ (12)

where: $U_{flow}^{max}$ = maximum flow velocity, Δt = time step, and Δx = spatial resolution, was maintained below the stability threshold (typically < 1 for explicit schemes). Given the grid resolution and expected maximum flow velocities within Polyfytos Lake (~0.5 m/s), typical grid cell size 200 m and time step 1 min (60 s), the resulting Courant number remained within the recommended safe range (CFL ≤ 0.15) throughout the domain.

C) Boundary conditions: Meteorological forcing was provided by the WRF model, which supplied spatially and temporally varying fields of wind speed, atmospheric pressure, air temperature, relative humidity, and solar radiation. These were mapped onto the Delft3D computational domain and formatted to meet model input specifications. Wind forcing was incorporated into the hydrodynamic model as shear stress acting on the free water surface [4]. The wind shear stresses in the x- and y-directions (τ_sx, τ_sy) were computed using:

${\tau }_{sx}={\rho }_a\cdot C_d\cdot U_{10m}^2\cos \left(\theta \right),{\tau }_{sy}={\rho }_a\cdot C_d\cdot U_{10m}^2\cos \left(\theta \right)$ (13)

where: ρ_α = air density (1.2 kg/m³), U_10m = wind speed at 10 m height and θ = wind direction angle relative to the model grid axes, both of them taken from the WRF model.

These stresses were converted to body forces $F_x^{wind}$ and $F_y^{wind}$ acting on the surface layer through:

$F_x^{wind}=\frac{1}{{\rho }_o}\cdot {\tau }_{sx},F_y^{wind}=\frac{1}{{\rho }_o}\cdot {\tau }_{sy}$ (14)

where: ρ_ο = density of water (1026.72 kg/m³). The drag coefficient C_d was calculated following the formulation by Smith and Banke [25]:

$C_d\left(U_{10m}\right)=\{\begin{array}{ll}{C}_d^A,\hfill & U_{10m}\le U_{10m}^A\hfill \\ C_d^A+\left(C_d^B-C_d^A\right)\frac{U_{10m}-U_{10m}^A}{U_{10m}^B-U_{10m}^A},\hfill & U_{10m}^A\le U_{10m}\le U_{10m}^B\hfill \\ C_d^B+\left(C_d^C-C_d^B\right)\frac{U_{10m}-U_{10m}^B}{U_{10m}^C-U_{10m}^B},\hfill & U_{10m}^B\le U_{10m}\le U_{10m}^C\hfill \\ C_d^C,\hfill & U_{10m}^C\le U_{10m}\hfill \end{array}$ (15)

where: $C_d^i$ = wind drag coefficients at the windspeed $U_{10m}^i$ ($i=A,B,C$ ) and are user defined. The values of these values are shown in Table 2.

Table 2. Meteorological stations adjacent to Lake Polyfytos.

User defined	$U_{10m}^i$ (m/s)	$C_d^i$
A	0	0.00063
B	10	0.00429
C	20	0.00495

The bottom friction in the Delft3D-FLOW model was parameterized using the Chezy formulation, a commonly used empirical approach to represent bed roughness effects in open channel and lake hydrodynamics. The Chezy coefficient C was set uniformly across the domain to:

• C = 65 for the U-direction, and

• C = 65 for the V-direction,

reflecting moderate bottom roughness conditions typical for a reservoir with mixed sediment and bedrock composition. The Chezy formulation expresses bed shear stress τ_b as:

${\tau }_b={\rho }_o\cdot g\cdot \frac{U_{d,avg}^2}{C^2}$ (16)

where: U_d,avg = depth-averaged velocity.

Lateral (wall) roughness was set to free-slip (wall roughness = 0), assuming minimal resistance from the lateral boundaries such as vertical walls or embankments, which is reasonable for the scale and resolution of the Polyfytos Lake model.

D) Processes: To resolve turbulent mixing processes in the hydrodynamic simulation, the model employed the standard k-ε turbulence closure scheme, which is suitable for simulating stratified flow conditions in large, deep-water bodies such as reservoirs and lakes. The following background (minimum) eddy viscosity and diffusivity values were prescribed to account for sub-grid mixing and ensure numerical stability in weakly turbulent regions:

• Horizontal eddy viscosity: 0.02 m²/s,

• horizontal diffusivity: 0.002 m²/s,

• vertical eddy viscosity: 2 × 10⁻⁵  m²/s,

• vertical diffusivity: 2 × 10⁻⁴  m²/s.

These values are within the typical range for large-scale inland water bodies and support the model’s ability to capture vertical stratification and mixing while minimizing numerical diffusion [4].

E) Heat flux model: Surface heat exchange was simulated using Delft3D-FLOW’s Ocean heat flux model, which accounts for various components of the surface energy budget including solar radiation, longwave radiation, sensible heat, and latent heat fluxes. To parameterize shortwave radiation attenuation, a Secchi depth of 5 meters was prescribed. This value represents the average water clarity in Polyfytos Lake and controls the vertical distribution of solar energy within the water column. The latent heat flux due to evaporation was calculated using the Dalton formula, with a Dalton number of 0.0013. Similarly, convective (sensible) heat exchange was governed by the Stanton number, also set to 0.0013 [4]. These coefficients represent typical values for moderately turbulent atmospheric conditions over large water bodies. Meteorological inputs required by the heat flux model, including: air temperature, relative humidity, and cloud coverage, were dynamically provided by the WRF atmospheric model, ensuring consistency between hydrodynamic and atmospheric boundary conditions.

F) Discharge Inflow and Outflow: To model the hydrodynamic behavior of Lake Polyfytos, both inflow and outflow discharge data were treated as critical boundary conditions in the Delft3D simulation. These discharges are primarily influenced by the operation of the Polyfytos Dam and hydroelectric plant, as well as the natural flow of the Aliakmon River. The inflow and outflow data were obtained from the work of Frysali et al. [16]. During the simulation period from 01 July to 08 July 2018, irrigation demand in the Thessaloniki area peaked at 30 m³/s, while precipitation was minimal—less than 2 mm—urther impacting the lake’s water balance.

Figure 8. Discharge inflow and outflow rates for Lake Polyfytos from 01 July to 07 July 2018.

Figure 8 displays the inflow and outflow discharge rates for Lake Polyfytos over this period, with time expressed in minutes on the x-axis. The inflow (blue line) remains relatively steady, fluctuating slightly between 1.5 and 2 m³/s, suggesting a consistent upstream water supply. In contrast, the outflow (orange line) shows a clear downward trend, starting at approximately 18 m³/s and gradually decreasing, with notable short-term increases observed around 7680 and 8640 minutes. This discrepancy between inflow and outflow suggests a net water loss from the lake during the simulated period, which likely impacts the lake’s water level and overall hydrodynamic behavior. These discharge values were incorporated into the Delft3D model as boundary conditions at the lake’s inflow and outflow grid points. Specifically, the inflow was applied at the bottom-left boundary of the model domain (see Figure 5), and the outflow at the top-right. This setup reflects the main flow direction observed in Lake Polyfytos.

3.2. Python Integration for Simulation Setup

Within the Delft3D interface, users can configure spatially varying wind and pressure conditions to support hydrodynamic and meteorological simulations. This interface serves as a comprehensive configuration tool, allowing users to define physical, numerical, and operational parameters relevant to a wide range of environmental modeling studies.

On the left panel and center of the interface (Figure 9), the integration with WRF is clearly visible. The Wind tab within the Physical parameters allows users to select wind and pressure fields that either vary spatially or remain uniform. These fields are likely sourced from WRF output, providing the flexibility to incorporate real-world atmospheric data into the simulation. Moving to the right panel, the Additional parameters section offers a detailed view of the interaction with Python code. Here, WRF-derived variables, such as #x_wind.amn# and #air_temperature.amf#, are mapped to internal keywords (e.g., Filwu, Filwt). This structure suggests a template-based configuration system, in which Python scripts likely parse these mappings and dynamically replace the tokens in simulation input files.

Figure 9. Delft3D interface for adding WRF parameters.

The Python code developed in this study is being released for the first time on GitHub under the name D3D_addin_parameters.py. This script is designed to significantly streamline the workflow for Delft3D users who wish to incorporate wind forcing and other atmospheric conditions derived from the WRF model. The code, provided in the appendix, supports a variety of meteorological variables, including:

• u10 represents the 10-meter U-wind component (in meters per second), while v10 represents the 10-meter V-wind component, both of which describe wind speed and direction at 10 meters above ground level,

• mslpressure is the mean sea level pressure, given in Pascals and converted to millibars for ease of interpretation. Similarly, surfpressure corresponds to the surface pressure, also converted from Pascals to millibars,

• tcc refers to total cloud cover, measured as a percentage, providing insight into the amount of cloud cover over the grid at a given time,

• t2m and d2m are the 2-meter air temperature and 2-meter dew point temperature, initially in Kelvin and converted to Celsius. These values help calculate other variables like relative humidity,

• msnswrf represents the shortwave radiation flux, measured in watts per square meter (W/m²), indicating the amount of solar radiation reaching the Earth’s surface,

• tp denotes total precipitation, initially in meters and converted to millimeters for consistency with common precipitation measurements,

• finally, mtpr gives the total precipitation rate in kilograms per square meter per second, converted to millimeters per hour for ease of use in simulation models.

Output Files Created:

1) air_temperature.amt: Contains air temperature (˚C).

2) relative_humidity.amr: Contains relative humidity (%).

3) air_pressure.amp: Contains air pressure (mbar).

4) x_wind.amu: Contains the x-component of wind speed (m/s).

5) y_wind.amv: Contains the y-component of wind speed (m/s).

6) cloudiness.amc: Contains cloudiness (%); note: this file was derived using a separate Python script (D3D_addin_CLD) and sourced from a .npz file format rather than WRF NetCDF outputs.

7) precipitation.ampr: Contains precipitation rate (mm/h).

8) sw_radiation_flux.ams: Contains shortwave radiation flux (W/m²).

3.3. Delft3D-PART Setup for Oil Spill Behavior

A) Description: Oil spill behavior in Polyfytos Lake was simulated using the PART module of Delft3D, which is specifically designed for particle tracking and allows simulation of pollutants—such as oil—only at the water surface and in Cartesian (x, y) coordinates.

B) Hydrodynamics: The modeling process began with the setup of the hydrodynamic conditions using the Delft3D-FLOW module. This involved defining the lake’s grid, bathymetry, and time-varying inflow and outflow discharge data. The Delft3D-FLOW model simulated both horizontal and vertical mixing, generating the current fields necessary for the subsequent oil transport simulations. Once the flow model was set up, the .com files were generated, which serve as fundamental inputs for the PART module.

C) Numerical parameters: A total of 400,000 particles were released into the model domain to represent the target substance(s), enabling detailed Lagrangian tracking of mass transport and spreading patterns over time. Vertical mixing of particles was modeled using the depth-averaged algebraic approach, which estimates vertical dispersion based on the resolved flow field and local turbulence levels. A scale factor of 0.1 was applied to modulate the strength of vertical dispersion, following Delft3D’s standard practice to tune sub-grid mixing. This approach is suitable for stratified conditions where vertical exchange is significantly weaker than horizontal dispersion and can be adjusted to reflect empirical or calibration data.

D) Continuous Releases: In the particle tracking simulation, a continuous oil release was modeled to represent a hypothetical surface spill scenario in Polyfytos Lake. The release was introduced at five predefined surface points: one located at the center of the lake, with two additional points positioned upstream and two downstream. These supplementary locations were spaced symmetrically at 3.4 km intervals from the central point, covering a total longitudinal span of approximately 13.6 km across the lake. Each release point was defined as a circular surface source with a radius of 10 meters, ensuring spatially distributed injection of particles to mimic realistic dispersion dynamics. At each location, the oil was released at a constant concentration of 20 kg/m³ and a steady flow rate of 0.1 m³/s, corresponding to a mass release rate of 2 kg/s per point. To adequately represent the distribution of the substance across all release points, 20% of the total 400,000 particles—equivalent to 80,000 particles-were allocated to each site. This configuration enabled the model to capture the influence of varying hydrodynamic conditions on the transport and fate of the oil across different regions of the lake.

The accuracy of a Lagrangian particle simulation is often expressed in terms of the smallest concentration that can be represented by a single particle within a computational cell [26]. This minimum resolvable concentration is calculated as the particle mass divided by the volume of the cell in which it resides:

$C_{min}=\frac{M_{total}}{N_{total}\cdot A_{cell}\cdot h_{layer}}$ (17)

where: C_min is the minimum concentration resolution (kg/m³), A_cell is the surface area of a computational (or zoom window) cell (m²), h_layer is the layer thickness (m) and M_total is the total mass released (kg), N_total = total number of released particles.

The above relationship can be used to estimate the required number of particles for a given simulation. In this study a hypothetical release of a total of 1,209,600 kg of substance. The smallest grid cell in the computational domain within the area of interest has a surface area of 40,000 m² with a minimum layer thickness of 0.1 m. It is required that the simulation should be able to resolve concentrations of at least 0.001 kg/m³. The minimum total number of particles required for this release would therefore be $N_{total}^{min}$ ≥ 1,209,600/(0.001 × 40,000 × 0.1) = 302,400 particles.

E) Process Parameters: To simulate subgrid-scale dispersion, Delft3D-PART [26] applies a Lagrangian random-walk scheme that combines deterministic advection with stochastic diffusion. The position of each particle is updated according to the following equation:

$x_i\left(t+\Delta t\right)=x_i\left(t\right)+U_{fvv}\cdot \Delta t+\eta \sqrt{6\cdot D\cdot \Delta t}$ (18)

where x_i(t) = position of the particle at time t, U_fvv = local flow velocity vector from Delft3D-FLOW, Δt = time step, D = turbulent dispersion coefficient and η = vector of independent random numbers drawn from the uniform interval [−1, 1].

This formulation reflects the combination of advective transport by resolved flow fields and stochastic displacement due to turbulent diffusion, with the random component scaled by $\sqrt{6\cdot D\cdot \Delta t}$ . The horizontal spreading of particles is governed by a time-dependent diffusion coefficient of the form:

$D=a\cdot t^b$ (19)

where α = 0.07 and b = 0.7 are calibration coefficients [26].

In the vertical settling/re-suspension equation from the Delft3D-PART manual [26] is:

${\upsilon }_s=c^n\left[A_0+A_1\sin \left(\frac{2\pi \left(t+\phi \right)}{T}\right)\right]$ (20)

where: υ_s is the settling (or re-suspension) velocity of particles, c = local particle concentration, n = exponent controlling concentration dependence, A₀ = baseline velocity component, A₁ = amplitude of the oscillating component, T = period of oscillation (e.g., daily or tidal cycle), φ = phase shift in time [same units as t], which adjusts the starting point of the sinusoidal variation.

In the present configuration, only the constant coefficient A₀ = 0.0001 m/s was applied, assuming a simplified vertical exchange driven by turbulence, without time-periodic modulation. Additionally, the exponent n was set to zero (n = 0), which implies cⁿ = 1, and settling velocity is independent of concentration.

F) Oil Model: The properties of the oil (e.g., density, viscosity) and the release timing were selected to reflect realistic environmental conditions, as illustrated in Figure 10.

In this simulation, we focus on the behavior of Heavy Fuel Oil (HFO), a dense, high-viscosity oil type known for its slow evaporation and strong adherence to surfaces. The configuration has been set with a detailed set of parameters to accurately represent how this oil would interact with the environment over time, specifically in a body of water like Polyfytos Lake. The following key parameters were derived from the works of Gonçalves et al. [27] and Tirado et al. [28].

Figure 10. Process parameters interface in Delft3D.

Key Parameters for Heavy Fuel Oil:

• Evaporation Rate: Set at 0.05 per day, the evaporation of Heavy Fuel Oil is slow, reflecting its nature as a less volatile substance. This is typical for heavy oils that do not quickly evaporate into the atmosphere.

• Stickiness Probability: With a value of 0.1, this indicates that the oil has a relatively low tendency to adhere to surfaces or particles, though it still maintains a moderate potential to stick, which could affect its spread across the water surface.

• Volatile Fraction: An unusually high value of 0.94 has been selected, indicating that 94% of the oil is capable of evaporating. It is acknowledged that the specified volatile fraction of 0.94 is uncharacteristically high for heavy fuel oil (HFO), which typically exhibits lower volatility. This parameter was initially selected to test upper-bound evaporation scenarios but will be re-evaluated and adjusted in future sensitivity analyses to reflect more realistic physical properties consistent with empirical HFO data. Such refinement is essential for ensuring the reliability of predicted oil weathering and transport dynamics.

• Emulsification Parameters: The emulsification behavior, crucial in understanding how the oil forms mixtures with water, has been set with a small emulsification rate of 2 × 10⁻⁶, and the maximum water content is set at 0.7. This suggests that the oil has a moderate potential to mix with water, forming emulsions, which could influence the oil’s movement and persistence on the water surface.

• Density: The density of the oil is 990 kg/m³, slightly lower than freshwater. This means that the oil will float on the water’s surface, spreading out in the form of a slick, but will also interact with the surface more easily than oils denser than water.

• Kinematic Viscosity: With a viscosity of 1500 centistokes (cSt), this oil is highly viscous, which means it will spread slowly across the surface of the water and resist rapid dispersion.

• Dispersion Settings: The model utilizes the Calculated dispersion option based on the Delvigne & Sweeney [29] approach. This allows the simulation to dynamically calculate how the oil particles disperse vertically due to turbulence in the water. This setting will adjust over time based on environmental conditions, influencing how the oil behaves in the water column.

Global Parameters:

• Minimum Slick Thickness: Set at 5 × 10⁻⁵ meters, this defines the threshold for detecting the slick, ensuring that only noticeable oil layers are tracked in the simulation.

• Deflection Angle (Coriolis Effect): The Coriolis effect is ignored (set to 0 degrees), a suitable choice for a small, inland water body like Polyfytos Lake, where this effect would be negligible.

4. Model Outputs and Results

The following analysis synthesizes model-generated outputs representing the spatial and temporal distribution of floating heavy fuel oil (measured in kg/m²) over a 7-day period. The data is derived from the Delft3D computational modeling, visualized through paired spatial maps for each key date range.

The model simulations for 01 July 2018 at 05:00 and 02 July 2018 at 11:00 reveal the initial release of oil into the reservoir. On the first day, the spill appears relatively contained, with low concentration levels ranging between approximately 1.36 to 2.73 kg/m², depicted in deep blue shades across most of the affected area. By the following day, a subtle expansion of the spill becomes evident—particularly toward the central region of the water body—signaling the onset of the spreading phase. This pair (Figure 11(a)) likely captures the very beginning of the spill event, marked by limited environmental influence. The prevailing wind conditions—light winds from the northwest to north-northwest at around 3 m/s—suggest minimal external forcing. As a result, the oil largely remains concentrated near its source, with little evidence of significant transport or widespread dispersion during this early stage.

Model outputs for 03 July 2018 at 17:00 (Figure 11(b)) indicate a sharp escalation in surface oil concentration. Distinct patches of green, yellow, and red hues emerge, representing localized zones with concentrations reaching up to approximately 13.64 kg/m². This surge in intensity appears to coincide with strong wind conditions—ranging from north-northeast to east, with speeds around 19 m/s—likely contributing to the abrupt increase in contamination. By 04 July 2018 at 23:00, while the most heavily contaminated zones begin to subside, the overall spatial extent of the spill expands considerably, particularly toward the northern

Figure 11. Temporal evolution of floating heavy fuel oil concentrations (kg/m²) in Polyfytos Lake from July 1 to July 7, 2018. Snapshots illustrate the spatial distribution and spread of the oil slick under varying wind conditions, with peak concentrations observed on July 3 during high wind events. Color gradients represent surface oil density, indicating both intensity and extent of dispersion.

and central regions of the water body. Although surface concentration levels decrease slightly, the geographic spread of the oil becomes more pronounced. This sequence reflects the continues release event, driven and intensified by high wind activity. The subsequent decrease in concentration, coupled with broader spatial dispersion, suggests the onset of natural degradation processes due to environmental dynamics.

Finally, the model outputs, from 06 July 2018 at 05:00 and 07 July 2018 at 11:00 (Figure 11(c)), highlight a continued westward expansion of the oil spill. On July 6th, the oil remains concentrated in the central to upper regions of the water body, with distinct red and yellow patches indicating localized areas of high surface density. By July 7th, the spill has both shifted and broadened, with its reach extending noticeably into the western sector. The dominant surface tones transition to cyan and green, corresponding to moderate concentration levels between 3 - 7 kg/m². High-density zones become more diffuse, suggesting active dispersion across a wider area. This final phase reflects ongoing horizontal spread, influenced by persistent winds from the north-northwest to northwest. The reduction in peak concentration levels, alongside an increase in the affected area, points to natural diffusion processes such as dispersion and evaporation.

5. Limitations and Uncertainties

This study has several limitations that merit discussion. First, the Delft3D and WRF models, while robust, rely on assumptions that may not fully capture small-scale turbulence, boundary interactions, or subsurface mixing processes. Parameterizations such as eddy viscosity and particle dispersion coefficients are generalized and not site-calibrated due to data constraints. Furthermore, uncertainties in input datasets—including bathymetric resolution, meteorological forecasts, and oil property specifications—can propagate through the simulation, affecting output reliability. Future work will focus on quantifying these uncertainties using ensemble simulations and sensitivity analyses.

Another limitation is the absence of direct comparison with real-world spill events or field trial data from Polyfytos Lake or similar inland water bodies. Incorporating such benchmarks would allow for more robust validation of modeled oil dispersion patterns and enhance the model’s operational relevance. This will be pursued in future studies, leveraging historical spill datasets where available or through collaboration with environmental monitoring agencies.

6. Conclusions

This study presents a comprehensive case of integrated hydrodynamic and atmospheric modeling applied to Polyfytos Lake in Greece, combining the Delft3D modeling suite with high-resolution meteorological data from the WRF model. Through the simulation of water circulation and the transport of a hypothetical oil spill, the research demonstrates the critical role of atmospheric forcing in determining pollutant dispersion in inland water bodies.

A key contribution of this work is the development and first-time deployment of a Python-based preprocessing pipeline that transforms WRF outputs into Delft3D-compatible meteorological input files. This codebase automates the conversion of critical atmospheric parameters—including wind components, temperature, humidity, pressure, cloudiness, radiation flux, and precipitation—into the required formats for spatially and temporally resolved modeling in Delft3D. This innovation significantly reduces setup time and technical complexity, offering a reusable and scalable workflow for future modeling efforts in different geographies or scenarios. Furthermore, a comprehensive set of steps within the PART (Particle Tracking) module is introduced.

Simulation results reveal that wind forcing is a dominant factor in shaping the spatial distribution and intensity of oil spills. Initially, the HFO remains concentrated near the release points, but as wind speed and direction change—especially with strong winds from the northeast—significant transport and spreading are observed. By July 3rd, oil concentration peaks at over 13 kg/m² in localized areas, emphasizing how meteorological dynamics accelerate surface dispersion. Subsequent model outputs show a gradual decline in surface concentration accompanied by horizontal diffusion toward the western and northern sections of the lake. This pattern suggests that natural degradation processes such as evaporation, dispersion, and emulsification begin to reduce the slick’s density while increasing its overall footprint.

Validation with local meteorological station data confirms that the WRF model accurately captures wind patterns and speed over the lake region, ensuring high-fidelity atmospheric forcing. Grid resolution tests indicate that while finer grids offer slightly more thermal precision (within 1.0% variation), a 200 m × 200 m resolution provides a robust trade-off between accuracy and computational efficiency, making it the optimal choice for large-scale operational modeling.

In summary, the case study illustrates how integrated hydrodynamic-atmospheric modeling, supported by tailored coding tools, can enhance decision-making for environmental monitoring, emergency response, and resource management. The methodology is not only applicable to Polyfytos Lake but is also transferable and reproducible for other inland water systems facing similar environmental challenges.

Acknowledgements

Part of this work, specifically the model output results, was previously included in a project deliverable submitted within the framework of the PERIVALLON project.

Author Contributions

V.P. conceived the study and wrote the abstract; V.P. and D.M. wrote the introduction and results; D.M. prepared the WRF atmospheric input data; C.A. prepared the figure illustrations; A.M. and K.V. reviewed the first draft; A.M., K.V., and I.G. supervised the methodology and conclusions sections; S.V. and I.K. acquired the funding. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

WRF output files used in this study are available from the corresponding author upon request. The scripts used to produce the appropriate files for the Delft3D input can be found at http://github.com/vaspapa79/D3D_usage.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Sharma, S., Gray, D.K., Read, J.S., O’Reilly, C.M., Schneider, P., Qudrat, A., et al. (2015) A Global Database of Lake Surface Temperatures Collected by in Situ and Satellite Methods from 1985-2009. Scientific Data, 2, Article No. 150008. https://doi.org/10.1038/sdata.2015.8
[2]	Li, X. and Tsigaris, P. (2024) The Global Value of Freshwater Lakes. Ecology Letters, 27, e14388. https://doi.org/10.1111/ele.14388
[3]	Hua, S., Jing, H., Qiu, G., Kuang, X., Andrews, C.B., Chen, X., et al. (2024) Long-Term Trends in Human-Induced Water Storage Changes for China Detected from GRACE Data. Journal of Environmental Management, 368, Article ID: 122253. https://doi.org/10.1016/j.jenvman.2024.122253
[4]	Hydraulics, D. (2024) User Manual of Delft3D-FLOW: Simulation and Multi-Dimensional Hydrodynamic Flows and Transport Phenomena, Including Sediments. WL/Delft Hydraulics. https://content.oss.deltares.nl/delft3d4/Delft3D-FLOW_User_Manual.pdf
[5]	Lesser, G.R., Roelvink, J.A., van Kester, J.A.T.M. and Stelling, G.S. (2004) Development and Validation of a Three-Dimensional Morphological Model. Coastal Engineering, 51, 883-915. https://doi.org/10.1016/j.coastaleng.2004.07.014
[6]	Curbani, F.E., Lacerda, K.C., Curbani, F., Barreto, F.T.C., Tadokoro, C.E. and Chacaltana, J.T.A. (2021) Numerical Study of Physical and Biogeochemical Processes Controlling Dissolved Oxygen in an Urbanized Subtropical Estuary: Vitória Island Estuarine System, Brazil. Environmental Modeling & Assessment, 27, 233-249. https://doi.org/10.1007/s10666-021-09787-1
[7]	Soulignac, F., Vinçon-Leite, B., Lemaire, B.J., Scarati Martins, J.R., Bonhomme, C., Dubois, P., et al. (2017) Performance Assessment of a 3D Hydrodynamic Model Using High Temporal Resolution Measurements in a Shallow Urban Lake. Environmental Modeling & Assessment, 22, 309-322. https://doi.org/10.1007/s10666-017-9548-4
[8]	Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Barker, D.M., Duda, M.G., Huang, X.Y., Wang, W. and Powers, J.G. (2008) A Description of the Advanced Research WRF Version 3 (No. NCAR/TN-475+STR). University Corporation for Atmospheric Research. http://dx.doi.org/10.5065/D68S4MVH
[9]	Anh, L.T., Anh, D.H., Yen Linh, M.T. and Thao, N.D. (2021) Investigation Typhoon Induced Storm Surge and High Wave in Vietnam Using Coupled Delft3d-Flow—WAVE Models Combined with Weather Research Forecast (WRF) Output Wind Field. Science & Technology Development Journal—Engineering and Technology, 4, 645-662. https://doi.org/10.32508/stdjet.v4i1.774
[10]	Guo, S., Zhu, D. and Chen, Y. (2023) Improvement and Evaluation of the Latest Version of WRF-Lake at a Deep Riverine Reservoir. Advances in Atmospheric Sciences, 40, 682-696. https://doi.org/10.1007/s00376-022-2180-5
[11]	Hassani, A., Santos, G.S., Schneider, P. and Castell, N. (2023) Interpolation, Satellite-Based Machine Learning, or Meteorological Simulation? A Comparison Analysis for Spatio-Temporal Mapping of Mesoscale Urban Air Temperature. Environmental Modeling & Assessment, 29, 291-306. https://doi.org/10.1007/s10666-023-09943-9
[12]	Li Ramírez, J.A., Zambrano Nájera, J.D.C. and Aristizábal Zuluaga, B.H. (2020) BVOC Emissions along the Eastern and Western Slopes of the Andes Central Range with Strong Altitudinal Gradient over a Wide Range of Andean Ecosystems: Model Estimation/Disaggregation with Biga. Environmental Modeling & Assessment, 25, 761-773. https://doi.org/10.1007/s10666-020-09698-7
[13]	Chen, X., Zhang, Y., Wang, Y., Yu, Q. and Ma, W. (2018) The Spatial-Scale Effect of an Atmospheric Environmental Impact Assessment in Regional Strategic Environmental Assessment (R-SEA). Environmental Modeling & Assessment, 23, 529-556. https://doi.org/10.1007/s10666-018-9590-x
[14]	Gabriel, M., Knightes, C., Cooter, E. and Dennis, R. (2014) The Impacts of Different Meteorology Data Sets on Nitrogen Fate and Transport in the SWAT Watershed Model. Environmental Modeling & Assessment, 19, 301-314. https://doi.org/10.1007/s10666-014-9400-z
[15]	Vatalis, K.I., Charalampides, G., Manoliadis, O., Akraios, C. and Asvesta, A. (2008) Water Quality Assessment of the Polyfyto Lake in Western Macedonia, Greece. Journal of Environmental Protection and Ecology, 9, 540-549. https://www.webofscience.com/wos/WOSCC/full-record/000260354500007
[16]	Frysali, D., Mallios, Z. and Theodossiou, N. (2023) Hydrologic Modeling of the Aliakmon River in Greece Using HEC-HMS and Open Data. Euro-Mediterranean Journal for Environmental Integration, 8, 539-555. https://doi.org/10.1007/s41207-023-00374-2
[17]	Papageorgiou, A., Papadakis, N. and Voutsa, D. (2015) Fate of Natural Organic Matter at a Full-Scale Drinking Water Treatment Plant in Greece. Environmental Science and Pollution Research, 23, 1841-1851. https://doi.org/10.1007/s11356-015-5433-3
[18]	Chortis, A.G., Ramandanis, I.G., Dadoulis, G.I., Melidis, L.N., Manos, G.C. and Katakalos, K.V. (2023). Structural Health Monitoring (SHM) of Polyfytos Bridge and Evaluation of Its Dynamic Properties. Proceedings of the 8th International Conference on Computational Methods in Structural Dynamics and Earthquake Engineering (COMPDYN 2015), Athens, 12-14 June 2023, 741-758. https://doi.org/10.7712/120123.10432.21194
[19]	Cantin, A., Beisner, B.E., Gunn, J.M., Prairie, Y.T. and Winter, J.G. (2011) Effects of Thermocline Deepening on Lake Plankton Communities. Canadian Journal of Fisheries and Aquatic Sciences, 68, 260-276. https://doi.org/10.1139/f10-138
[20]	Morrison, H., Curry, J.A. and Khvorostyanov, V.I. (2005) A New Double-Moment Microphysics Parameterization for Application in Cloud and Climate Models. Part I: Description. Journal of the Atmospheric Sciences, 62, 1665-1677. https://doi.org/10.1175/jas3446.1
[21]	Clough, S.A., Shephard, M.W., Mlawer, E.J., Delamere, J.S., Iacono, M.J., Cady-Pereira, K., et al. (2005) Atmospheric Radiative Transfer Modeling: A Summary of the AER Codes. Journal of Quantitative Spectroscopy and Radiative Transfer, 91, 233-244. https://doi.org/10.1016/j.jqsrt.2004.05.058
[22]	Betts, A.K. and Miller, M.J. (1986) A New Convective Adjustment Scheme. Part II: Single Column Tests Using GATE Wave, BOMEX, ATEX and Arctic Air-Mass Data Sets. Quarterly Journal of the Royal Meteorological Society, 112, 693-709. https://doi.org/10.1002/qj.49711247308
[23]	Mellor, G.L. and Yamada, T. (1982) Development of a Turbulence Closure Model for Geophysical Fluid Problems. Reviews of Geophysics, 20, 851-875. https://doi.org/10.1029/rg020i004p00851
[24]	Lagouvardos, K., Kotroni, V., Bezes, A., Koletsis, I., Kopania, T., Lykoudis, S., et al. (2017) The Automatic Weather Stations NOANN Network of the National Observatory of Athens: Operation and Database. Geoscience Data Journal, 4, 4-16. https://doi.org/10.1002/gdj3.44
[25]	Smith, S.D. and Banke, E.G. (1975) Variation of the Sea Surface Drag Coefficient with Wind Speed. Quarterly Journal of the Royal Meteorological Society, 101, 665-673. https://doi.org/10.1002/qj.49710142920
[26]	Deltares (2024) Delft3D-PART User Manual: Particle Tracking Module for Delft3D Suite. https://content.oss.deltares.nl/delft3d4/Delft3D-PART_User_Manual.pdf
[27]	Santos, R.G., Loh, W., Bannwart, A.C. and Trevisan, O.V. (2014) An Overview of Heavy Oil Properties and Its Recovery and Transportation Methods. Brazilian Journal of Chemical Engineering, 31, 571-590. https://doi.org/10.1590/0104-6632.20140313s00001853
[28]	Tirado, A., Félix, G., Trejo, F., Varfolomeev, M.A., Yuan, C., Nurgaliev, D.K., Sámano, V. and Ancheyta, J. (2023) Properties of Heavy and Extra-Heavy Crude Oils. In: Ancheyta, J., Varfolomeev, M. and Yuan, C., Eds., Catalytic In-Situ Up-Grading of Heavy and Extra-Heavy Crude Oils, Wiley, 1-38. https://doi.org/10.1002/9781119871507.ch1
[29]	Delvigne, G.A.L. and Sweeney, C.E. (1988) Natural Dispersion of Oil. Oil and Chemical Pollution, 4, 281-310. https://doi.org/10.1016/s0269-8579(88)80003-0