An Assessment of Computational Fluid Dynamics as a Tool to Aid the Design of the HCMR-Artificial-ReefsTM Diving Oasis in the Underwater Biotechnological Park of Crete

Since recreational diving activities have increased in recent decades, resulting in additional environmental pressure on the coastal zone, the deployment of artificial reefs as a conservation strategy to divert mass ecotourism from fragile natural reefs has been proposed and realized in many areas of the world. Twelve units of a patented naturoid artificial reef technology developed by the Hellenic Centre for Marine Research (HCMR) were deployed in 2015 in the Underwater Biotechnological Park of Crete (UBPC) in order to create an experimental diving oasis and investigate the potential of achieving this aim for the over-exploited coastal ecosystems of this part of the Eastern Mediterranean. Assessment of the degree of establishment of artificial reefs and their ability to mimic natural ecosystems is often monitored through biological surveys and sampling. The measurement of the chemical, physical, and hydrodynamic characteristics of the water mass surrounding artificial reefs is also essential to fully understand their comparison to natural reefs. In particular, the flow field around reefs has been shown to be one of the most important physical factors in determining suitable conditions for the establishment of a number of key species on reef habitats. However, the combination of biological establishment monitoring and realistic flow-field simulation using computational fluid dynamics as a tool to aid in the design improvement of already existing reef installations has not been fully investigated in previous work. They are often reported separately as either ecological or engineering studies. Therefore, this study examined a full-scale numerical simulation of the field flow around individual already installed naturoid reef shapes, and part of their present arrangement on the sea bottom of the UPBC combined with the field-testing of the functionality of the installed artificial reefs concerning fish species aggregation. The results show that the simulated flow characteristics around the HCMR diving oasis artificial reefs were in good general agreement with the results of former studies, both for flows around a single deployed unit and for flows around a cluster of more than one unit. The results also gave good indications of the performance of individual reef units concerning key desirable characteristics such as downstream shadowing and sediment/nutrient upwelling and resuspension. In particular, they confirmed extended low flow levels (less than 0.3 m/s) and in some cases double vortexes on the downstream side of reef units where observed colonization and habitation of some key fish species had taken place. They also showed how the present distribution of units could be optimized to perform better as an integrated reef cluster. The use of computational fluid dynamics, with field survey data, is therefore suggested as a useful design improvement tool for installed reef structures and their deployment arrangement for recreational diving oases that can aid the sustainable development of the coastal zone.

of natural interest in order to relieve manmade pressure on their local ecosystems. The investigation was carried out in the context of the five-year field-testing results from the biological survey and sampling of the deployed ARs in order to evaluate the combination of CFD and in situ AR establishment monitoring as a design improvement tool for the HCMR diving oasis. In particular, this research examined the additional benefit of carrying out a full-scale numerical simulation of the actual naturoid reef shapes, and part of their present arrangement on the sea bottom of the UBPC, along with the field-testing of the functionality of the UBPC artificial reefs concerning fish species aggregation. Fish aggregation around AR units is reported in these studies to be driven by field flow characteristics. The combination of the upstream upwelling flow with the formation of a wake zone downstream of an AR unit results in sediment/nutrient resuspension and the presence of eddies and vortexes. These flow fields provide shelter, feeding ground, spawning ground, and rest areas for fish assemblages around an AR [70]. Consequently, this research is presented both as an evaluation and future-design improvement tool while possible multidisciplinary uses and benefits of such an application are also discussed.

Study Area
In the framework of various research projects of the Institute of Marine Biology, Biotechnology and Aquaculture (IMBBC) of the Hellenic Centre for Marine Research (HCMR), the Underwater Biotechnological Park of Crete (UBPC) has been developed. The UBPC is a unique large-scale in situ research infrastructure that covers an area of 30,000 m 2 , lies about 2 km offshore the northern Cretan coast in the vicinity of the HCMR land premises in Crete, with a depth increasing from 18 to 22 m along the South-North direction (Figure 1). The objectives of the UBPC are the protection, conservation, and exploitation of marine biological resources and the development of innovative technologies in the fields of open sea invertebrate aquaculture and marine ecotourism, including the development and testing of artificial habitats, as well as long-term monitoring of the coastal environment. In May 2014, a bottom moored seafloor observatory was deployed that is comprised of various scientific instruments mounted on a metal truss construction ( Figure 2). The seafloor observatory is equipped with an Acoustic Doppler Profiler (ADP) (for the measurement of sea current velocity/direction in the water column, wave velocity, and height), a fluorometer (for the measurement of in-vivo chlorophyll-a, turbidity, and fluorescence) and a conductivity-temperature-density instrument (CTD, for the continuous measurement of sea temperature, salinity, density, sound velocity, and dissolved oxygen). Alongside the seafloor observatory, a vertical temperature logger array is deployed that has six sensors equally spaced from 7 m down to 17 m depth.
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 25 survey and sampling of the deployed ARs in order to evaluate the combination of CFD and in situ AR establishment monitoring as a design improvement tool for the HCMR diving oasis. In particular, this research examined the additional benefit of carrying out a full-scale numerical simulation of the actual naturoid reef shapes, and part of their present arrangement on the sea bottom of the UBPC, along with the field-testing of the functionality of the UBPC artificial reefs concerning fish species aggregation. Fish aggregation around AR units is reported in these studies to be driven by field flow characteristics. The combination of the upstream upwelling flow with the formation of a wake zone downstream of an AR unit results in sediment/nutrient resuspension and the presence of eddies and vortexes. These flow fields provide shelter, feeding ground, spawning ground, and rest areas for fish assemblages around an AR [70]. Consequently, this research is presented both as an evaluation and future-design improvement tool while possible multidisciplinary uses and benefits of such an application are also discussed.

Study Area
In the framework of various research projects of the Institute of Marine Biology, Biotechnology and Aquaculture (IMBBC) of the Hellenic Centre for Marine Research (HCMR), the Underwater Biotechnological Park of Crete (UBPC) has been developed. The UBPC is a unique large-scale in situ research infrastructure that covers an area of 30,000 m 2 , lies about 2 km offshore the northern Cretan coast in the vicinity of the HCMR land premises in Crete, with a depth increasing from 18 to 22 m along the South-North direction (Figure 1). The objectives of the UBPC are the protection, conservation, and exploitation of marine biological resources and the development of innovative technologies in the fields of open sea invertebrate aquaculture and marine ecotourism, including the development and testing of artificial habitats, as well as long-term monitoring of the coastal environment. In May 2014, a bottom moored seafloor observatory was deployed that is comprised of various scientific instruments mounted on a metal truss construction ( Figure 2). The seafloor observatory is equipped with an Acoustic Doppler Profiler (ADP) (for the measurement of sea current velocity/direction in the water column, wave velocity, and height), a fluorometer (for the measurement of in-vivo chlorophyll-a, turbidity, and fluorescence) and a conductivity-temperaturedensity instrument (CTD, for the continuous measurement of sea temperature, salinity, density, sound velocity, and dissolved oxygen). Alongside the seafloor observatory, a vertical temperature logger array is deployed that has six sensors equally spaced from 7 m down to 17 m depth. The aforementioned instruments that comprise the seafloor observatory assist the AR infrastructure since the continuous logging of environmental data is vital when it comes to the monitoring of marine organisms-how they cope with the inter-annual seasonal changes and determining the impact of climatic extreme conditions, minima and maxima-especially for a study area in the northern Cretan coast that is prone to Lessepsian migration [71,72].

Constructing and Deploying the HCMR Artificial Reefs
The construction of the AR units took place at the land premises of the HCMR in Gournes, Crete (Figure 3a), which are approximately 2 km from their deployment area (Figure 3b). The ARs were made of concrete using plastic modules as the construction's shaping frame. The selection of concrete had to do with the desired resemblance to natural rocky reef habitats (Figures 4a and 4b), not only in an aesthetic pattern but also in the functionality of their surface. The concrete's rough texture provides "microhabitats" and larger refugia more prone to biofouling than the ones of a smooth material surface and eventually lead to a more naturoid design [73]. Finally, a pair of stainless steel pipes were used to make the transport mounts for every unit. These pipes run parallel to each other through the base of every unit and provide an ergonomic design for safe transport and sea deployment operations. Lift belts for cranes are easily mounted on four slots for every AR during transport and are easily removed by divers after the deployment on the seafloor. Furthermore, every AR unit had a 100 mm vertical hole through its center axis in case an additional mounting pipe was going to be inserted in the seafloor to provide more stability to counteract scouring effects or trawling activity. However, for this study, it was decided not to insert additional mounting pipes in the seafloor, in order to observe the operational behavior of the AR and how they responded both to natural (bad weather conditions, scouring effects) and manmade (bottom trawling, fishing nets, anchors) interactions. The aforementioned instruments that comprise the seafloor observatory assist the AR infrastructure since the continuous logging of environmental data is vital when it comes to the monitoring of marine organisms-how they cope with the inter-annual seasonal changes and determining the impact of climatic extreme conditions, minima and maxima-especially for a study area in the northern Cretan coast that is prone to Lessepsian migration [71,72].

Constructing and Deploying the HCMR Artificial Reefs
The construction of the AR units took place at the land premises of the HCMR in Gournes, Crete (Figure 3a), which are approximately 2 km from their deployment area (Figure 3b). The ARs were made of concrete using plastic modules as the construction's shaping frame. The selection of concrete had to do with the desired resemblance to natural rocky reef habitats (Figure 4a,b), not only in an aesthetic pattern but also in the functionality of their surface. The concrete's rough texture provides "microhabitats" and larger refugia more prone to biofouling than the ones of a smooth material surface and eventually lead to a more naturoid design [73]. Finally, a pair of stainless steel pipes were used to make the transport mounts for every unit. These pipes run parallel to each other through the base of every unit and provide an ergonomic design for safe transport and sea deployment operations. Lift belts for cranes are easily mounted on four slots for every AR during transport and are easily removed by divers after the deployment on the seafloor. Furthermore, every AR unit had a 100 mm vertical hole through its center axis in case an additional mounting pipe was going to be inserted in the seafloor to provide more stability to counteract scouring effects or trawling activity. However, for this study, it was decided not to insert additional mounting pipes in the seafloor, in order to observe the operational behavior of the AR and how they responded both to natural (bad weather conditions, scouring effects) and manmade (bottom trawling, fishing nets, anchors) interactions.
The constructed ARs were separated into two different groups concerning their structural geometric morphology. The first group (Type A in Figure 4c) provides big vertical crevices, most of which are blind (do not continue all the way through the construct), thus allowing bigger distinct external and internal surfaces for every AR unit. The second group (Type B in Figure 4c) provides big horizontal crevices-both blind and through-giving this group of ARs up to four distinct vertical levels of concrete. AR units of both groups weigh between 1500 and 2500 kg, have a height between 1.4 and 2.5 m, and have base diameters between 1.2 and 2.3 m. In addition, a small AR unit, with code B1 in Figure 4c, was constructed for use as a checkpoint through the diving path of the UBPC, and weighed roughly 400 kg. AR unit had a 100 mm vertical hole through its center axis in case an additional mounting pipe was going to be inserted in the seafloor to provide more stability to counteract scouring effects or trawling activity. However, for this study, it was decided not to insert additional mounting pipes in the seafloor, in order to observe the operational behavior of the AR and how they responded both to natural (bad weather conditions, scouring effects) and manmade (bottom trawling, fishing nets, anchors) interactions.   The constructed ARs were separated into two different groups concerning their structural geometric morphology. The first group (Type A in Figure 4c) provides big vertical crevices, most of which are blind (do not continue all the way through the construct), thus allowing bigger distinct external and internal surfaces for every AR unit. The second group (Type B in Figure 4c) provides big horizontal crevices-both blind and through-giving this group of ARs up to four distinct vertical levels of concrete. AR units of both groups weigh between 1500 and 2500 kg, have a height between 1.4 and 2.5 m, and have base diameters between 1.2 and 2.3 m. In addition, a small AR unit, with code B1 in Figure 4c, was constructed for use as a checkpoint through the diving path of the UBPC, and weighed roughly 400 kg. Prior to the deployment operation of the ARs was the in situ survey conducted by the members of the Scientific Diving Team of the HCMR. The local characteristics of the UBPC were recorded through many transects and the morphology and seagrass coverage of the seafloor were taken into account for the careful selection of the exact places that the ARs were going to be placed. These transects were carried out on a seasonal basis and were realized using underwater recordings, both photos and video. Visual material was obtained either by SCUBA operations or by the deployment of stand-alone cameras. For the post-sampling analysis and for divers' convenience, checkpoint marks (coded floats) were deployed on the seafloor to create a grid both for spatial reference and for navigational aid. The final AR placement selection was with respect to the local habitat that is mainly muddy sand, covered by meadows of the green macroalgae Caulerpa prolifera (Lamouroux, 1809) interspersed with small patches of Posidonia oceanica (Delile, 1813).  Prior to the deployment operation of the ARs was the in situ survey conducted by the members of the Scientific Diving Team of the HCMR. The local characteristics of the UBPC were recorded through many transects and the morphology and seagrass coverage of the seafloor were taken into account for the careful selection of the exact places that the ARs were going to be placed. These transects were carried out on a seasonal basis and were realized using underwater recordings, both photos and video. Visual material was obtained either by SCUBA operations or by the deployment of stand-alone cameras. For the post-sampling analysis and for divers' convenience, checkpoint marks (coded floats) were deployed on the seafloor to create a grid both for spatial reference and for navigational aid. The final AR placement selection was with respect to the local habitat that is mainly muddy sand, covered by meadows of the green macroalgae Caulerpa prolifera (Lamouroux, 1809) interspersed with small patches of Posidonia oceanica (Delile, 1813).

Creating the Computational Grid for 3-D Demonstration and Computational Fluid Dynamics (CFD) Analysis
Once the AR units were constructed, photos of them were taken, as in Figure 5a, and by using the Autodesk 123D CAD and modeling tools software suite, a 3-D mesh for every unit was created, as in Figure 5b. Once the geometry was extracted, the 3-D mesh provided the realization of various demonstrational alternatives. The HCMR-AR TM concept was presented to various audiences to raise awareness of sustainable and marine-organisms-oriented coastal zone development through many educational and citizen science initiatives of the HCMR. Not only were recorded videos of the AR development timeline presented, but also 3-D printed models of the actual units were utilized for better understanding of the purpose of these submerged structures that attract fish and divers' interest. Once the AR units were constructed, photos of them were taken, as in Figure 5a, and by using the Autodesk 123D CAD and modeling tools software suite, a 3-D mesh for every unit was created, as in Figure 5b. Once the geometry was extracted, the 3-D mesh provided the realization of various demonstrational alternatives. The HCMR-AR TM concept was presented to various audiences to raise awareness of sustainable and marine-organisms-oriented coastal zone development through many educational and citizen science initiatives of the HCMR. Not only were recorded videos of the AR development timeline presented, but also 3-D printed models of the actual units were utilized for better understanding of the purpose of these submerged structures that attract fish and divers' interest. The extracted mesh file of every AR unit was imported in ANSYS Workbench software. Therein, the ICEM package was used to set up and prepare the computational mesh of different AR units for CFD analysis. Following the ANSYS tutorials and guidelines [74], a water tunnel was designed with dimensions of 48 m × 24 m × 12 m (length, width, and height respectively) to study one single AR unit. The under-investigation AR unit's center was placed approximately 12 m from the inlet wall, for the horizontal distance to be 5 times the height of the AR-object. In a similar way, the height of the water tunnel was defined as 12 m and the width of the tunnel was defined as 24 m. In these terms, the AR unit's mesh was at least 5 times its height away from any wall surface of the tunnel. For the lee side of the AR unit, the distance to the outlet was greater than 15 times its height (~36 m) ( Figure  6a). For the study of an AR cluster, the dimensions of the simulation water tunnel were 60 m × 30 m × 15 m (length, width, and height respectively) in order to keep the aforementioned protocols concerning distancing-AR height ratios (Figure 6b). The unstructured mesh that was generated was then imported into ANSYS Fluent for CFD analysis. The k-ε model is widely used for practical engineering modeling of turbulent flows due to its robustness, economy, and reasonable accuracy for a wide range of flows [60,62,65]. The k-ε model assumes that the flow is fully turbulent and that the effects of molecular viscosity are negligible. This pressure-based k-ε RNG method was selected. For the solution methods, the Green-Gauss cell-based was selected for the estimation of the gradients. The second-order upwind was selected for momentum, turbulent kinetic energy, turbulent dissipation rate, and energy. For the pressure-velocity coupling, the SIMPLEC algorithm scheme was selected. The temperature was set to 293 K (which was the average value measured by the instruments of the UBPC); the initial gauge pressure was set to 300,000 Pa to simulate the almost 20 The extracted mesh file of every AR unit was imported in ANSYS Workbench software. Therein, the ICEM package was used to set up and prepare the computational mesh of different AR units for CFD analysis. Following the ANSYS tutorials and guidelines [74], a water tunnel was designed with dimensions of 48 m × 24 m × 12 m (length, width, and height respectively) to study one single AR unit. The under-investigation AR unit's center was placed approximately 12 m from the inlet wall, for the horizontal distance to be 5 times the height of the AR-object. In a similar way, the height of the water tunnel was defined as 12 m and the width of the tunnel was defined as 24 m. In these terms, the AR unit's mesh was at least 5 times its height away from any wall surface of the tunnel. For the lee side of the AR unit, the distance to the outlet was greater than 15 times its height (~36 m) (Figure 6a). For the study of an AR cluster, the dimensions of the simulation water tunnel were 60 m × 30 m × 15 m (length, width, and height respectively) in order to keep the aforementioned protocols concerning distancing-AR height ratios (Figure 6b). The unstructured mesh that was generated was then imported into ANSYS Fluent for CFD analysis. The k-ε model is widely used for practical engineering modeling of turbulent flows due to its robustness, economy, and reasonable accuracy for a wide range of flows [60,62,65]. The k-ε model assumes that the flow is fully turbulent and that the effects of molecular viscosity are negligible. This pressure-based k-ε RNG method was selected. For the solution methods, the Green-Gauss cell-based was selected for the estimation of the gradients. The second-order upwind was selected for momentum, turbulent kinetic energy, turbulent dissipation rate, and energy. For the pressure-velocity coupling, the SIMPLEC algorithm scheme was selected. The temperature was set to 293 K (which was the average value measured by the instruments of the UBPC); the initial gauge pressure was set to 300,000 Pa to simulate the almost 20 m depth of the actual conditions; and the top surface was defined as a moving wall with zero initial velocity and no-slip shear conditions. With respect to the data recorded from the ADP of the seafloor observatory of the UBPC, the maximum velocity for the first four meters above the seabed did not exceed 60 cm/s, thus the input velocities for the CFD analysis were defined as 0.5 m/s and 1 m/s. The lower value (0.5 m/s) was representative of the actual higher current velocities recorded in the UBPC and the higher value (1 m/s) was selected as a safety factor indication of a really extreme flow speed scenario for our study area.

Monitoring of the Deployed AR
Immediately following the deployment of the artificial reefs (ARs) in the UBPC, monitoring

Monitoring of the Deployed AR
Immediately following the deployment of the artificial reefs (ARs) in the UBPC, monitoring dives took place to record the structural integrity and the interaction of every unit with the natural environment, both in technical and ecological terms. The first individuals of damselfish-Chromis chromis (Linnaeus, 1758)-settled in the reefs during the first week after the deployment (Figure 7a). In the third week after the AR deployment, some large individuals of groupers-family Serranidae, Epinephelus costae (Steindachner, 1878)-occupied large cavities of the units (Figure 7b). Six weeks after the deployment, the first individuals of the family Sparidae and especially white seabream-Diplodus sargus (Linnaeus, 1758)-had already settled in the vertical crevices of various ARs (Figure 7c,d). In addition, in the seventh week after the AR deployment, the first individuals of dusky groupers-Epinephelus marginatus (Lowe, 1834)-were attracted and already settled in the AR crevices ( Figure 7e). Additionally, in less than two months after the deployment of the ARs the first pelagic hunters-greater amberjacks, Seriola dumerilii (Risso, 1810)-were encountered in large schools of a few hundreds of individuals (Figure 7f). Videos and additional information can be found on the website of the UBPC: http://ubpcrete.hcmr.gr/index.php/timeline. Epinephelus costae (Steindachner, 1878)-occupied large cavities of the units (Figure 7b). Six weeks after the deployment, the first individuals of the family Sparidae and especially white seabream-Diplodus sargus (Linnaeus, 1758)-had already settled in the vertical crevices of various ARs (Figures  7c and 7d). In addition, in the seventh week after the AR deployment, the first individuals of dusky groupers-Epinephelus marginatus (Lowe, 1834)-were attracted and already settled in the AR crevices ( Figure 7e). Additionally, in less than two months after the deployment of the ARs the first pelagic hunters-greater amberjacks, Seriola dumerilii (Risso, 1810)-were encountered in large schools of a few hundreds of individuals (Figure 7f). Videos and additional information can be found on the website of the UBPC: http://ubpcrete.hcmr.gr/index.php/timeline. A major concern for the whole project was the fact that the study area of the UBPC is in open sea and easily accessible by artisanal fishermen and spearfishing divers. Although a serious effort and a campaign took place to inform the local fishermen about the research activities and the (a) (c) A major concern for the whole project was the fact that the study area of the UBPC is in open sea and easily accessible by artisanal fishermen and spearfishing divers. Although a serious effort and a campaign took place to inform the local fishermen about the research activities and the demonstrational purposes of the deployed AR, the scientific divers of the UBPC found that three units were turned over and discovered remaining fishing lines and nets on various ARs (Figure 8a), proving that these units were pulled to rotate and fall down. As an advantage to the design of the ARs comes the fact that even when some ARs were pulled over onto the seafloor, they still retained their functionality as an attractive habitat for various species and their crevices still hosted organisms that used to dwell there before the aforementioned manmade interventions. Following further testing of the turned-over ARs, concerning their long-term functionality, it was decided not to rotate them back to their initial vertical position, in order to witness any possible decline in marine organisms' aggregation. For the time being, no radical changes occurred and the "fallen" ARs still retain their functionality, as the dorsal fins of a "hiding" dusky grouper shows in Figure 8b, still occupying its habitat. From Figure 8a, it is also noticeable that there was a similarity of the AR texture and appearance with a natural rocky reef, as the external surface of the AR was completely covered by a dense community of seaweeds. Two years after the deployment of the AR, a baseline survey of sessile benthos was carried out [75]. A non-destructive assessment method was applied by means of photoquadrat sampling by the HCMR Scientific Diving Team. Twenty-five quadrats (0.25 × 0.25 m) were randomly obtained from five AR units (five replicates per reef) and were analyzed. The study indicated a high colonization potential for sessile benthic taxa, when considering the biotic coverage and number of taxa developed on the artificial surfaces, especially in the absence of natural hard substrata in the broader area. Figure A1 in the Appendix provides information on substrate cover percentages from the aforementioned study.
Sustainability 2020, 12, x FOR PEER REVIEW 10 of 25 of the turned-over ARs, concerning their long-term functionality, it was decided not to rotate them back to their initial vertical position, in order to witness any possible decline in marine organisms' aggregation. For the time being, no radical changes occurred and the "fallen" ARs still retain their functionality, as the dorsal fins of a "hiding" dusky grouper shows in Figure 8b, still occupying its habitat. From Figure 8a, it is also noticeable that there was a similarity of the AR texture and appearance with a natural rocky reef, as the external surface of the AR was completely covered by a dense community of seaweeds. Two years after the deployment of the AR, a baseline survey of sessile benthos was carried out [75]. A non-destructive assessment method was applied by means of photoquadrat sampling by the HCMR Scientific Diving Team. Twenty-five quadrats (0.25 × 0.25 m) were randomly obtained from five AR units (five replicates per reef) and were analyzed. The study indicated a high colonization potential for sessile benthic taxa, when considering the biotic coverage and number of taxa developed on the artificial surfaces, especially in the absence of natural hard substrata in the broader area. Figure A1 in the Appendix provides information on substrate cover percentages from the aforementioned study. During the last five years that followed the deployment of the ARs, a wide variety of marine organisms attractive to recreational divers were encountered and recorded in the artificial oasis of the UBPC, including octopuses, moray eels, Mediterranean cardinal fish, goldblotch groupers, dusky groupers, triggerfish, greater amberjacks, tuna, European barracuda, sea bream, dentex, loggerhead turtles, etc.

Computational Fluid Dynamics (CFD) Simulation of a Single AR Unit
Selected artificial reefs (ARs) were simulated and the computational results are presented. For a (a) (b) Figure 8. Four years after the deployment of the ARs. (a) A unit covered with biofouling organisms fully integrated with the natural "seascape" and fishing nets and lines. (b) As a result of fishing activity this unit was overturned, but still in its crevices a dusky grouper is hiding and finds shelter, as it is revealed from the presence of its dorsal fins.
During the last five years that followed the deployment of the ARs, a wide variety of marine organisms attractive to recreational divers were encountered and recorded in the artificial oasis of the UBPC, including octopuses, moray eels, Mediterranean cardinal fish, goldblotch groupers, dusky groupers, triggerfish, greater amberjacks, tuna, European barracuda, sea bream, dentex, loggerhead turtles, etc.

Computational Fluid Dynamics (CFD) Simulation of a Single AR Unit
Selected artificial reefs (ARs) were simulated and the computational results are presented. For a better demonstration of the 3-D flow characteristics, the results are presented in a horizontal and a vertical plane. The horizontal plane was set at half the height of the AR and the vertical plane was set as a cross-section of the middle of the AR, parallel to the flow input direction. The ARs of this study were not symmetrical since they were not designed to resemble a defined geometric shape, but to mimic natural rocky reefs. Consequently, the complexity of the generated mesh and of the actual flow made the simulations very demanding concerning the available computational resources.
In Figure 9a, the vertical cross-section of the flow around an AR unit (with code A5-refer to Figure 4c) is presented. The input velocity was 0.5 m/s and in front of the AR, the flow velocity decreased because of the AR surface resistance. Upwelling flows occured and the maximum velocities were obtained in the highest area of the unit. In the upstream region, there were small areas with various increasing velocities. Downstream of the AR, a wake region was defined and a large vortex formed. The wake region expanded more than 3 m in the horizontal axis, thus creating the desired shade effect for fish aggregation there, since the flow velocity was less than 0.3-0.35 m/s. In Figure 9a, the vertical cross-section of the flow around an AR unit (with code A5-refer to Figure 4c) is presented. The input velocity was 0.5 m/s and in front of the AR, the flow velocity decreased because of the AR surface resistance. Upwelling flows occured and the maximum velocities were obtained in the highest area of the unit. In the upstream region, there were small areas with various increasing velocities. Downstream of the AR, a wake region was defined and a large vortex formed. The wake region expanded more than 3 m in the horizontal axis, thus creating the desired shade effect for fish aggregation there, since the flow velocity was less than 0.3-0.35 m/s. In Figure 9b, the horizontal plane of the same velocity (0.5 m/s) flow simulation is presented. The plane was located half the height of the AR unit above the bottom surface, to be less affected by the boundary interactions between the flow and the bottom surface. In the upstream region, the velocity again decreased in front of the AR unit due to its resistance, and flow velocity increased in the flanks of the AR unit, where maximum values were obtained. A complex vortex was generated in the downstream region and the desired conditions of flow velocity (less than 0.3 m/s) were realized In Figure 9b, the horizontal plane of the same velocity (0.5 m/s) flow simulation is presented. The plane was located half the height of the AR unit above the bottom surface, to be less affected by the boundary interactions between the flow and the bottom surface. In the upstream region, the velocity again decreased in front of the AR unit due to its resistance, and flow velocity increased in the flanks of the AR unit, where maximum values were obtained. A complex vortex was generated in the downstream region and the desired conditions of flow velocity (less than 0.3 m/s) were realized in an area that exceeded 3 m.
In Figure 10, the flow simulation of the same AR unit is presented, with flow velocity input equal to 1 m/s. Similar to the 0.5 m/s analysis, the same principles and results governed this numerical experiment. In the vertical cross-section, the same flow patterns were visible. The inlet velocity declined in the upstream region in front of the AR unit and the maximum values were obtained in a narrow area upstream of the AR unit. Upwelling fields were witnessed again, but the relative intensity of the flow velocity increase appeared to be lesser when compared with the lower inlet flow velocity, hence the denser red vectors above the AR unit in Figure 9a. In the downstream area, a wake region that expanded more than 3 m horizontally was created. A vortex was generated in the wake area downstream of the AR unit and, again, this area provided the desired shelter for fish to rest and aggregate (Figure 10a). Concerning the horizontal cross-section of the 1 m/s experiment (Figure 10b), there was a decrease in the upstream region velocity and the maximum values were obtained in the flanks of the AR unit. The wake area was also visible and the vortex and shade-effect area expanded up to 3 m.
Sustainability 2020, 12, x FOR PEER REVIEW 12 of 25 experiment. In the vertical cross-section, the same flow patterns were visible. The inlet velocity declined in the upstream region in front of the AR unit and the maximum values were obtained in a narrow area upstream of the AR unit. Upwelling fields were witnessed again, but the relative intensity of the flow velocity increase appeared to be lesser when compared with the lower inlet flow velocity, hence the denser red vectors above the AR unit in Figure 9a. In the downstream area, a wake region that expanded more than 3 m horizontally was created. A vortex was generated in the wake area downstream of the AR unit and, again, this area provided the desired shelter for fish to rest and aggregate (Figure 10a). Concerning the horizontal cross-section of the 1 m/s experiment (Figure 10b), there was a decrease in the upstream region velocity and the maximum values were obtained in the flanks of the AR unit. The wake area was also visible and the vortex and shade-effect area expanded up to 3 m. In Figures A2 and A3 of the Appendix, the respective flow simulations with an inlet flow velocity of 0.5 m/s and 1 m/s of the B3 AR unit (refer to Figure 4c) are presented. These numerical simulations resulted in similar flow field characteristics to the ones presented here for the A5 AR unit. In Figures A2 and A3 of the Appendix, the respective flow simulations with an inlet flow velocity of 0.5 m/s and 1 m/s of the B3 AR unit (refer to Figure 4c) are presented. These numerical simulations resulted in similar flow field characteristics to the ones presented here for the A5 AR unit.

Computational Fluid Dynamics (CFD) Simulation of an Artificial Reef (AR) cluster
This study included the flow simulation around a cluster of ARs. This numerical analysis was not carried out to exactly simulate the flow conditions in the UBPC but aimed to provide an insight into the interaction between a part of the arrangement of the submerged structures of the diving oasis and the flow in terms of an AR unit being in the proximity of another. Due to the arrangement configuration, mutual affected flow characteristics occur and may impact any desired effects. For example, in Figure 11a seven cross-section planes vertical to the flow direction were presented for the flow-interactions between ARs to be visible. A 0.2 m height horizontal plane was also presented to assist in the interpretation of the flow fields. The vertical cross-cut planes were equally spaced every 2.5 m and the shadow effect in the downstream region behind a unit could be seen. Additionally, in the cross-cuts shown in Figure 11a(ii,iv), the decrease in flow velocity upstream of the ARs is visible. The increase of velocity in the flanks and upper region of every AR unit can also be seen in the respective cross-cut planes with yellowish contours (Figure 11a(i,iii,v)).

Computational Fluid Dynamics (CFD) Simulation of an Artificial Reef (AR) cluster
This study included the flow simulation around a cluster of ARs. This numerical analysis was not carried out to exactly simulate the flow conditions in the UBPC but aimed to provide an insight into the interaction between a part of the arrangement of the submerged structures of the diving oasis and the flow in terms of an AR unit being in the proximity of another. Due to the arrangement configuration, mutual affected flow characteristics occur and may impact any desired effects. For example, in Figure 11a seven cross-section planes vertical to the flow direction were presented for the flow-interactions between ARs to be visible. A 0.2 m height horizontal plane was also presented to assist in the interpretation of the flow fields. The vertical cross-cut planes were equally spaced every 2.5 m and the shadow effect in the downstream region behind a unit could be seen. Additionally, in the cross-cuts shown in Figure 11a (ii) and (iv), the decrease in flow velocity upstream of the ARs is visible. The increase of velocity in the flanks and upper region of every AR unit can also be seen in the respective cross-cut planes with yellowish contours (Figure 11a i,iii,v).   In Figure A4 of the Appendix, the respective CFD simulation with an inlet flow velocity of 1 m/s for the same AR cluster shown in this section is presented. This numerical simulation resulted in similar flow field characteristics to the ones presented here.

Discussion
Many local authorities and stakeholders have already expressed their interest in developing a recreational facility such as the HCMR diving oasis and the benefits of such a project are numerous. From an ecological point of view, the installation of oases in environmentally degraded areas has been proved to relieve areas of natural interest (MPAs, etc.) from the ever-growing recreational diving pressure and protect and support of the local biodiversity in the neighboring areas of oases, as well as provide more effective supervision and control of diving activities if they are facilitated in such defined areas. From a socio-economic point of view, in terms of integrated coastal zone management, a diving oasis can offer benefits that enable engaged co-existence between all end-users of the coastal zone involved, from divers and fishermen to local businesses and tourism infrastructure. A further benefit is the opportunity for marine education and training activities, as well as the raising of local public awareness about various environmental topics. However, the HCMR diving oasis is a system at a demonstration level requiring continuous evaluation in terms of design and functionality performance.
As part of this evaluation, this study was an assessment of a numerical simulation tool that employs computational fluid dynamics to calculate the flow around the naturoid artificial reefs (ARs) already deployed as the HCMR experimental diving oasis in the UBPC. This was combined with a brief timeline of the progress of fish establishment that are attractive to divers. The assessment was therefore in terms of the simulations being used as a design and functionality improvement tool for already deployed ARs and the results from this AR functionality assessment are in the context of the experimental diving oasis having already been realized in the field in 2015. Further background that should be remembered is that the diving oasis area was of no recreational diving interest before installation of the ARs, and within two months became attractive to target species that may interest a diver with these species very quickly aggregating and finding permanent shelter in the AR infrastructure.
The numerical analysis presented in this paper was in good agreement with previous studies and the results of the flow simulation proved the AR potential to create the desired shade effect in the downstream side of an AR unit. This shade effect provided sufficient space for turbulent flow conditions with velocities of up to 0.35 m/s, which has been reported to be the preferred range for fish to rest and seek shelter [65]. Moreover, the numerical analysis was also in good agreement with previous studies in that it showed that the main flow characteristics of the HCMR-AR TM were the upwelling field and the reverse flow in the vortex-dominated wake area [61]. Additionally observed was that the height of an AR unit may increase the slow-flow zone's scale and the back vortex flow's intensity behind it [61,65], something that in our study was shown when comparing the flow fields of two different AR units (Figure 9a compared to Figures A2a and 10a compared to Figure A3a) under the same inlet flow velocity. Additionally, when comparing the same AR units/cluster with different inlet velocities, the flow field characteristics did not change qualitatively, but the length of the back-vortex flow decreased with increasing incoming flow velocity [61]. Moreover, the desired upwelling flow characteristic for an AR unit was presented in this study. An upwelling field is defined where the flow velocity component of the z-axis (height) is equal to or greater than 10 percent of inlet flow velocity [60]. Upwelling is reported by numerous studies to increase the bottom and surface water sediment and nutrients transport [59,60,63]. By this resuspension, upwelling enhances biological productivity, which feeds fish, thereby creating a habitat richer in nutrients and oxygen diffusion. The consequent thriving of algae and plankton are reported by Liu et al. [61] to be the driving force for the reefs not only to provide shelter from predators but also an abundant food source for aggregating fish. In our simulations, the upwelling fields were presented above the AR units with the red vectors (Figures 9a, 10a, A2a and A3a) facing towards the upper water layers after the flow separation in the top surface of every unit. The area affected by the upwelling seemed to increase with the increase in the height and/or the incident flow area of the respective AR (e.g., Figure 9a compared to Figure A2a), a result that is in good agreement with the ones of former studies. Jiang et al. showed that the insertion of a guide plate to disturb the flow around a cubic AR had a positive overall function in altering the flow fields around an AR [62]. In the future, some of our new units could be designed to adopt the functionality of such a guide plate while still retaining their naturoid design concept. This could possibly have a positive impact for AR units that are going to be deployed in areas of low current speeds, where poor sediment-nutrients resuspension is expected. The understanding and the deliberate alteration of the upwelling field can only be calculated and realized through the utilization of CFD and the parameterization of key input variables (AR dimensions, current speed/orientation, etc.). When considering this and the following discussion, it should be noted that, facing the same computational problems as former studies, the numerical simulations here simplified the seabed topography and represented it as a flat surface with no indents or formations [62,65].
The simulation for the flow characteristics around an AR cluster presented in this study, rather than a single unit, was different in some aspects when compared with other studies. One main difference resulted from the diving oasis concept of our study when compared with other AR sets, concerning the reef unit deployment arrangement. Since the main purpose for the deployment of other AR sets was not the creation of a recreational diving facility, but to enhance the local marine livestock through fish aggregation around the AR sets, different practices were followed. Firstly, the AR units that were deployed in other areas were designed not to be diver-attractive, so they have defined geometric shapes (e.g., cubes) since they do not attempt to mimic natural rocky reefs. Secondly, since the AR positioning on the seafloor of other studies was not part of a diving path, different AR units may be either have been piled up one on another or AR clusters may even have had a distancing of up to 100 m or more [12]. For our diving oasis concept, the AR units had to be close enough to each other for the next unit to be visited to be easily discernible for a diver, in order to follow a diving path. This prerequisite led to a spacing of maximum of 10 m between reefs, thus a distance of almost five times the height of a typical unit that is 2 m in height was possible. According to our results, with increasing inlet flow velocity, the length of the back-vortex flow decreased. This suggests that the spacing should be reduced in order to increase the intensity of the desired shadow effect in cases like Figure 11b and also maintain its characteristics as much as possible in higher flow velocities. This is something that is in good agreement with the findings of previous studies which show that the impact of two AR units is strongest when the reef units spacing is between three and four times their height [61]. Worth mentioning here is the fact that various other studies define the ARs according to the legislation of their local authorities, e.g., Korean regulations where an AR set should have a facility volume greater than 800 m 3 . Typical AR units used in these studies are cubic-shaped with 2 m edges (volume of 8 m 3 ), resulting in more than 100 units to create an AR set [63]. Such a number of deployed AR units was outside the scope of our diving oasis concept, since a number of 15-20 single units per diving trail seemed sufficient.
Concerning the flow characteristics of the individual AR units, the numerical simulation presented here proved to be a valuable tool for their functional assessment and their future design details and confirmed the AR units' functionality. The desired upwelling and shadowing effects in the upstream and the downstream area, respectively, of every single unit were witnessed through the CFD simulations. However, the results also suggested that the shape of the AR units could be altered in order to increase the area of the vortex dominated wake region. This is important when the direction of the flow is perpendicular to the gap between two AR units, where the interaction of the two flow fields was almost negligible as shown in our study (Figure 11a) and in previous work [61,65]. These previous studies with cubic reef units have shown that when spacing increases to two reef heights, the flow field presents a separation trend.
The numerical simulation results here thus provide vital information for the operation of future AR units in the diving oasis. The shape, openings, and orientation of every future AR unit are going to be defined and pre-simulated in order to enhance the desired effects and characteristics seen in these results. Furthermore, the colonization of key species that was encouraged by the designed flow characteristics of the ARs is going to be witnessed with the installation of a live network camera in the UBPC that will provide live-streaming images under any weather/sea-state condition [76]. Thus, the imperceptible to the human eye flow domains, which affect fish behavior, will be studied with the visualization of the flow through numerical simulations. Using the CFD tool tested here, the possible impacts of AR geometry and construction variations (including crevices, protuberances and their size) on key fish species may thus be validated concerning different flow conditions and various unit deployment arrangements.
Furthermore, the results here highlight the necessity of the CFD analysis as a tool for the design improvement of the experimental diving oasis concerning the very nature of the AR units. Since the units here were constructed to mimic the sparse rocky reef natural habitat of the northern coast of Crete, the assessment and prediction of their flow fields were only achievable through numerical modeling. Any new units will have a large number of openings and crevices of different diameters that need to be recalculated, not only for the technical integrity of the new ARs [59], but also for increasing habitat availability for the functionality of their carrying capacity of crevice-dwelling organisms as artificial habitats [77]. In the present study, the geometry of the ARs were simplified and the horizontal crevices of the Type B units were simulated as blind, since the scope of this study was the investigation of the characteristics of the flow field around a whole AR unit and clusters of units. There is therefore further scope for the development of this application of the CFD tool for a more detailed micro-analysis of the flow fields. Furthermore, the development of an experimental water tunnel for flow simulations is scheduled for the near future. This water tunnel, assisted by 3-D printing technology, will implement particle image velocimetry (PIV) techniques that, accompanied by the CFD simulation techniques developed here, will provide a more accurate and holistic approach to the R&D projects of the flow fields around the submerged structures. Moreover, the further development of CFD could support numerous other research projects of the UBPC (e.g., studying the impact of the flow field on the growth rate/shape of sessile invertebrates of experimental aquaculture [78]). Finally, it is worth mentioning the necessity for scientific environmental monitoring, beyond the in situ scientific diving surveys of the colonization of the ARs. These provide essential data on the bio-chemical and physical characteristics of the surrounding sea water that inform us about the ecological footprint of any manmade deployed structure. The environmental monitoring at the UBPC utilizes in situ scientific loggers (e.g., CTDs, ADPs, temperature logger arrays) to give data on such parameters as temperature, salinity, and current speed and direction. For our study location, available remote sensing datasets were compared with in situ autonomous temperature logger datasets resulting in highly correlated match-ups [79], thus providing confidence in using remote sensing data as a proxy for ecological monitoring and also as a source for former years' local environmental and climatic information. Further work on the integration of these environmental variables should assist future numerical simulations of the conditions around the AR and the diving oasis to be as accurate and analytical as possible, comprising and modeling all the available input parameters.

Conclusions
The AR design of the HCMR diving oasis aimed at (a) diver-attractiveness, (b) mimicking the natural environment, and (c) fish aggregation, thus defined geometry shapes were excluded while the parameterization in our studies included various other inputs. For example, different overall AR unit shapes/crevices, different arrangement/distancing, and different orientation to the local prevailing current conditions needed to be taken into consideration. These could only be fully realized with the aid of CFD numerical simulations of the flow fields, in order to deliver locally tailored AR units and AR arrangements that combined the aforementioned proposed innovative characteristics into a sustainable recreational diving infrastructure. Overall, the numerical simulations were in good qualitative agreement with the ones of former studies, since the expected field flow characteristics were observed, concerning the regimes both in the upstream and the downstream region of AR units.
Specifically, the utilization of CFD, when combined with in situ data (environmental, field surveys, and visual census), was shown to be a useful tool in providing key information on the performance of the naturoid-type HCMR-ARs TM in the diving oasis field setting. CFD analysis proved to be a valuable tool not only for optimizing the design of one AR unit but also in determining the placement and the relative position and distancing between the units of an AR cluster and thus the overall design of a diving oasis. For future designs and deployments, the local geomorphological and environmental characteristics and prevailing sea current conditions should also be further incorporated. Once they are fully combined with the CFD modeling, this will allow the incorporation of any AR constructions into a diving oasis in the most locally tailored manner possible and under the scope of developing eco-friendly manmade diving infrastructures that can improve the sustainability of the coastal zone.