**Table 1.** List of investigated types of flour, together with their identification code (ID) and typology.

### *2.2. Gamma Spectrometry Analysis and Radiological Hazard Effects Assessment*

For the gamma spectrometry investigations, samples were packed in a 20 mL polyethylene plastic vial to reach a geometry homogeneity around the detector; then, the respective net weights were measured and recorded. After this preparation step, they were counted for 70,000 s, and spectra were analyzed in order to obtain the activity concentration of 137Cs and 40K, through the evaluation of their γ-lines at 661.66 keV and 1460.81 keV, respectively. A positively-biased Ortec HPGe detector (GEM) (FWHM of 1.85 keV; peak-to-Compton ratio of 64:1; relative efficiency of 40% at the 1.33 MeV 60Co γ-line), placed inside lead wells to shield the background radiation environment, was employed for the analysis. A photo of the experimental setup is reported in Figure 2.

**Figure 2.** HPGe gamma spectrometry setup.

A multipeak Marinelli geometry gamma source (AK-5901) of 1 L capacity, covering the energy range 59.54–1836 keV, customized to reproduce the exact geometries of samples in a water-equivalent epoxy resin matrix, was employed for efficiency and energy calibrations. The ANGLE 4 code was employed for the efficiency transfer factor calculations to the 20 mL vial sample holder geometry [33]. For data acquisition and analysis, Gamma Vision (Ortec) software was used [34].

The activity concentration (*C*) of each detected radionuclide was calculated according to the following [35,36]:

$$\mathcal{C} = \frac{N\_E}{\varepsilon\_E t \gamma\_d M} \tag{1}$$

where *NE* indicates the net area of the radioisotope photopeak; *εE* and *γd* are the efficiency and yield of the photopeak, respectively; *M* is the mass sample (g); and *t* is the lifetime (s).

The resulting uncertainty of each measurement, with coverage factor *k* = 2, was calculated, taking into account the uncertainty of the counting estimation, of the calibration source, of the efficiency calibration, of the background subtraction, and of the γ-branching ratio [37].

The quality of the gamma spectrometry experimental results was certified by the Italian Accreditation Body (ACCREDIA) [38].

This implies the continued verification (with annual periodicity) of the maintenance of the gamma spectrometry method performance characteristics.

The repeatability of the results, in particular, is verified over time with the double test method. A certified reference material (also containing the radionuclides 40K and 137Cs) is analyzed twice, with the active concentration of the radionuclide of interest defined as *x*1 (first measurement) and *x*2 (second measurement). The probability level *p* = 0.95 is considered. The following formula is applied:

$$|x\_1 - x\_2| \le \sqrt{2} \cdot s\_{r'} t^r$$

where *t* is the Student variable, and *sr* is the standard deviation of repeatability obtained in the validation phase.

To evaluate radiological hazard effects, the effective dose for flour ingestion was calculated [39]:

$$D\_{\rm ing} \left( \frac{\rm Sv}{\rm y} \right) = h\_{\rm ing, K-40} \times f\_{\rm ing, K-40} \tag{2}$$

where *hing,K-40* is the coefficient of the effective dose for the insertion unit for ingestion of 40K (Sv Bq−1), and *Jing,K-40* is the intake of 40K (Bq year−1) [40]. The latter value is obtained by multiplying the yearly flour consumption for the mean activity concentration of the investigated radionuclide experimentally measured.

### *2.3. ICP-MS and ICP-OES Analysis*

For the analysis, approximately 0.5–1.0 g of sample, together with 5 mL of ultrapure (for trace analysis) HNO3 (67–69%) and 1 mL of hydrogen peroxide (H2O2), as well as 5 mL of ultrapure water, were directly introduced into a 100 mL TFM vessel. Acid digestion was performed using a Milestone microwave system, Ethos 1, in four steps: 5 min at 1000 W and 85 ◦C, 10 min at 1000 W and 145 ◦C, 4 min at 1000 W and 180 ◦C, 15 min at 1000 W, and cooling [41]. The mixture was filtered and filled up to 100 mL with distilled H2O.

ICP-MS calibration solutions for the analytes measured at 1, 5, 10, 100, and 200 ppb for As, Be, Cd, Cu, Pb, and Zn, as well as at 0.2, 0.5, 1.0, 2.0, and 3.0 ppb for Hg, were prepared by opportune dilutions of two different reference materials (o2si Multi-element aqueous RM 100 mL 20 μg mL−1, and o2si Single-element Mercury Hg aqueous RM 100 mL 1007 μg mL−1) in 0.5 % (*v*/*v*) HNO3 and 0.5% (*v*/*v*) HCl [42]. Internal standardization was applied, with an internal standard at 20 ppb, added on-line via a T-piece before the nebulizer. For the ICP-MS measurements a Thermo Scientific iCAP Qc ICP-MS was used [43] (see Figure 3).

The sample introduction system consisted of a cooled (3 ◦C) Peltier, a baffled cyclonic spray chamber, a PFA nebulizer, and a quartz torch with a 2.5 mm i.d. removable quartz injector. The instrument operated in a single collision cell mode, with kinetic energy discrimination (KED), using pure He as collision gas. All samples were presented for analysis using a Cetac ASX-520. The iCAP Qc ICP-MS operated in a single KED mode using 1550 W forward power, 0.98 L min−1 nebulizer gas, 0.8 L min−1 auxiliary gas, 14.0 L min−1 cool gas flow, 4.5 mL min−1 He collision cell gas, 45 s each for sample uptake/wash time, optimized dwell times per analyte (0.1 s, except 0.5 s for As), one point per peak, and three repeats per sample.

**Figure 3.** ICP-MS experimental setup.

For the ICP-OES measurements, a Perkin Elmer Optima 2100 DV was employed [44] (see Figure 4).

**Figure 4.** ICP-OES experimental setup.

The operating conditions were 1450 W RF power, 15 L min−1 plasma gas flow rate, 0.2 L min−1 auxiliary gas flow rate, 0.8 L min−1 nebulizer gas flow rate, 1.5 L min−1 sample flow rate, 50 s read delay, CCD detector, and nitrogen purge gas, with 589.592 nm (Na), 285.213 nm (Mg), 317.933 nm (Ca) wavelengths.

### *2.4. Multivariate Statistical Analysis*

The Microcal Origin software for Windows (Northampton, MA, United States, ver. 2018) was used for all statistical calculations [45]. A statistical analysis was conducted, starting on a multivariate matrix, where variables were the average concentrations of detected elements (P, Mg, Ca, Na, Fe, Zn, Mn, Cu, Ni), and the cases were the eleven analyzed flour samples.

With the aim of unambiguously correlating flour samples to their botanical origin, according to their elemental concentrations, an exploratory method (principal components analysis, PCA) was performed. It ensures the reduction of the data dimensionality, whereas the combinations of variables identified by the PCs provide the greatest contribution to sample variability. Before the elaboration, the appropriateness of the data set was

preliminarily checked by the Bartlett test (sphericity test), as already reported in a previous work [46].

### **3. Results and Discussion**

*3.1. Radioactivity Analysis*

The activity concentration of 40K in the investigated flour samples, as calculated by Equation (1), is reported in Table 2 as mean value ± standard deviation for each type of flour.

The highest value refers to chestnut flour. This may reasonably reflect the different composition of the two investigated flours, reported to have 847 mg and 171 mg of potassium per 100 g of sample [47,48], respectively.

Regarding 137Cs, its activity concentration in all analyzed samples, reported in Table 2, turned out to be lower than the minimum detectable activity value, thus excluding an anthropogenic radioactive contamination of the investigated samples.

### *3.2. The Annual Effective Dose for Food Ingestion*

Table 3 reports the estimation of the annual effective dose due to the ingestion of investigated samples, *Ding*, as calculated by Equation (2), for the age category higher than 17 years.

**Table 3.** Annual effective dose due to the ingestion of investigated samples, *Ding*, for the age category higher than 17 years, together with its percentage with respect to worldwide average natural dose to humans (2.4 mSv y−1).

Reported values take into account the average yearly consumption of flour in Italy, assumed to be of about 85 kg, taking into consideration its use to produce bread, pizza, and baked goods [49,50], making the hypothesis that this need was satisfied from a single type of flour. For each type of investigated flour, Table 3 also shows the percentage values of the aforementioned doses with respect to the total (external + internal) natural radioactivity value for humans, expressed in terms of effective dose, equal to 2.4 mSv y−1 [51]. Noteworthily, the calculated annual effective doses for all investigated flour types were under the threshold level of 1 mSv y−1 imposed by the legislation [51] (see Figure 5).

**Figure 5.** *Ding* for all investigated samples, compared to the threshold level imposed by legislation.

Therefore, no risk to humans due to ionizing radiations can be associated with the flour ingestion.

### *3.3. Heavy Metal Analysis*

Table 4 reports Cd and Pb mean concentrations (mg kg−1) for each type of investigated flour, as obtained through ICP-MS analysis.

**Table 4.** Cd and Pb mean concentrations in investigated flour samples.

As can be seen from the table, the experimental values are lower than the contamination thresholds established by the commission regulations [52]. Therefore, Cd and Pb cannot be considered as pollutants, and are generally expected not to cause severe effects on human health.

### *3.4. Elemental and Chemometric Analysis*

Table 5 reports the elemental composition of the analyzed flour samples.

**Table 5.** Average mineral content, determined by ICP-OES and ICP-MS, for each investigated flour type.

In particular, macro- (P, Mg, Ca) and micro-inorganic minerals (Na, Fe, Zn, Mn, Cu, Ni) determined by ICP-OES and ICP-MS were tabulated; their average concentrations were in good agreement with values reported in literature [53].

As reported [48], the presence of inorganic minerals in flour samples is strictly correlated with the plant or fruit of origin. In order to unambiguously establish this correlation, PCA was conducted, considering nine variables and eleven cases. Preliminarily, the measure of sampling adequacy and Bartlett's test of sphericity were carried out to verify the suitability of the data for factor analysis. Test results suggest that the correlation matrix was factored and appropriate for PCA. By using the Kaiser criterion, two principal components (PC1 and PC2) with eigenvalues exceeding one were extracted, representing the 40.57% and 23.72% of total variance, respectively.

PCA results reveal the presence of four different groups, reported in the bi-plot of Figure 6.

**Figure 6.** Results of the PCA (bi-plot) for all analyzed samples.

The groups appear to be comparable to each other in terms of macro- and microelement concentration. Specifically, four groups were recognized: group I (containing samples S10 and S11, green points) and IV (containing samples S1 and S8, blue points) refer to flours that turned out to be positively correlated, as far as PC1 is concerned, with Ca, Mg, P, Cu, Zn, Fe, Ni, and Mn, and negatively correlated with Na. Regarding PC2, flours in group I exhibited a positive correlation with values of Cu, P, Mg, and Ca contents, and a negative correlation with other variables. The opposite occurred for flours in group IV. Continuing on, group II (containing samples S4, S5, and S9, purple points) and III (containing samples S2, S3, S6, and S7, black points) refers to flours that turned out to be negatively correlated, as far as PC1 is concerned, with Ca, Mg, P, Cu, Zn, Fe, Ni, and Mn, and positively correlated with Na. Regarding PC2, flours in group II exhibited a positive correlation with the values of the Cu, P, Mg, and Ca contents, and a negative correlation with other variables. The opposite occurs for flours in group III.

The differences in the average concentrations of P, Mg, Ca, Na, Fe, Zn, Mn, Cu, and Ni revealed among the identified groups strictly depend on the compositional features of the plant or fruit from which the analyzed flours derive. In the light of this, the obtained results could potentially be helpful to the producers to warn against food fraud, ensure flour quality, and inform the customers about sensitive parameters to be used potentially in the traceability process, in order to build up several fingerprints able to ensure the healthy effects of these products.

### **4. Conclusions**

In this article, radiation levels, metals contamination, and inorganic mineral content in flours were measured; this is extremely important, due to their high consumption by the population for nutritional use.

The specific activity of the main natural radionuclide present in flour, 40K, was measured using HpGe gamma spectrometry, with the aim of estimating the health risk for the age category higher than 17 years by the assessment of the effective dose due to flour ingestion. The coefficient of the effective dose for ingestion reported by Italian Legislative Decree 101/20 was employed. The estimation of the annual effective dose was found to be in the range 0.04–0.23 mSv y−1, about 1.6–9.6% of the total natural radioactivity value for humans (2.4 mSv y−1) and under the threshold level (1 mSv y−1), thus excluding the risk of ionizing radiation effects on humans. The 137Cs activity concentration was lower than the minimum detectable activity value for all analyzed samples, thus excluding its presence as pollutant.

As far as the heavy metal (Cd and Pb) concentration is concerned, it was lower than the contamination threshold value, thus excluding their presence as contaminants.

The PCA multivariate statistical analysis, performed starting from the inorganic minerals content in the investigated samples, showed how it is possible to discriminate among them, with a strict dependance on the compositional features of the plant or fruit from which the analyzed flours derive, thus representing a powerful tool to guarantee the product authenticity to the consumers.

**Author Contributions:** Conceptualization, F.C. and V.V.; methodology, F.C. and V.C.; validation, D.M.; formal analysis, G.A., M.D., A.B., S.M. and M.M.; investigation, F.C., G.A., G.P. and V.V.; resources, F.C., V.C. and D.M.; data curation, F.C.; writing—original draft preparation, F.C.; supervision, D.M., V.C. and V.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

### *Article* **Assessment of Natural Radioactivity and Radiological Risks in River Sediments from Calabria (Southern Italy)**

**Francesco Caridi 1,\*, Marcella Di Bella 2,3, Giuseppe Sabatino 1, Giovanna Belmusto 4, Maria Rita Fede 5, Davide Romano 1,2, Francesco Italiano 2 and Antonio Francesco Mottese 1**

**Abstract:** This study was developed to carry out a comprehensive radiological assessment of natural radioactivity for river sediment samples from Calabria, southern Italy, and to define a baseline background for the area on a radiation map. In the studied area, elevated levels of natural radionuclides are expected, due to the outcropping acidic intrusive and metamorphic rocks from which the radioactive elements derive. To identify and quantify the natural radioisotopes, ninety river sediment samples from nine selected coastal sampling points (ten samples for each point) were collected as representative of the Ionian and the Tyrrhenian coastline of Calabria. The samples were analyzed using a gamma ray spectrometer equipped with a high-purity germanium (HPGe) detector. The values of mean activity concentrations of 226Ra, 232Th and 40K measured for the studied samples are (21.3 ± 6.3) Bq kg−1, (30.3 ± 4.5) Bq kg−1 and (849 ± 79) Bq kg−1, respectively. The calculated radiological hazard indices showed average values of 63 nGy h−1 (absorbed dose rate), 0.078 mSv y−1 (effective dose outdoors), 0.111 mSv y−1 (effective dose indoors), 63 Bq kg−1 (radium equivalent), 0.35 (Hex), 0.41 (Hin), 0.50 (activity concentration index) and 458 μSv y−1 (Annual Gonadal Equivalent Dose, AGED). In order to delineate the spatial distribution of natural radionuclides on the radiological map and to identify the areas with low, medium and high radioactivity values, the Surfer 10 software was employed. Finally, the multivariate statistical analysis was performed to deduce the interdependency and any existing relationships between the radiological indices and the concentrations of the radionuclides. The results of this study, also compared with values of other locations of the Italian Peninsula characterized by similar local geological conditions, can be used as a baseline for future investigations about radioactivity background in the investigated area.

**Keywords:** natural radioactivity; river sediments; gamma spectroscopy; radiological risks; background radioactivity; multivariate statistics

### **1. Introduction**

Radioactivity naturally occurs in the environment. Aside from the cosmic radiation which bombards the earth from the outside of the atmosphere, the predominant natural source of human exposure to radiation doses is due to the terrestrial radiation which includes radionuclides of 238U, 232Th, 235U series and 40K [1].

It is established that the abundance of these primordial radionuclides depends on the local environment's geology [2,3]. Moreover, natural and anthropogenic events such

**Citation:** Caridi, F.; Di Bella, M.; Sabatino, G.; Belmusto, G.; Fede, M.R.; Romano, D.; Italiano, F.; Mottese, A.F. Assessment of Natural Radioactivity and Radiological Risks in River Sediments from Calabria (Southern Italy). *Appl. Sci.* **2021**, *11*, 1729. https://doi.org/10.3390/ app11041729

Academic Editor: Richard Kouzes

Received: 20 January 2021 Accepted: 11 February 2021 Published: 15 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

as flooding, radionuclide uptake by plants and the use of fertilizers, usually provoke an alteration in the distribution of the abundance levels of the radionuclides in the environment [4–6].

Studies of natural radiation are needed to establish reference levels, especially in areas where the hazard of radioactive exposure may be higher due to the lithological features of the geological context.

In this framework, central and southern Calabria (Southern Italy) can be considered as "risky" areas, since granitoids and metamorphic rocks that represent sources of radionuclides widely outcrop along the Serre and the Aspromonte Massifs. Moreover, recent studies highlighted the presence of heavy minerals, such as monazite and zircons, in the Calabrian rocks [7], which might contain high amounts of radioactive U and Th.

In this study, fluvial deposits from nine different selected sites from the Ionian and Tyrrhenian coastline of Calabria were analyzed. River sediments, consisting of mineral particles with different size and chemical composition, are considered reliable long-term indicators of river pollution by radionuclides, because water pollution components are deposited in the sediments [8]. Major sources of natural radionuclides in sediments have different possible origins, like weathering and recycling by erosion of terrestrial minerals and rocks (igneous or metamorphic) containing 40K and radionuclides of uranium and thorium radioactive series, and also rainfall and other depositional phenomena such as landslides and precipitations [9]. The main aim of this paper is to identify and quantify natural gamma-emitting radionuclides and to evaluate the radiological hazard effects [10]. For this purpose, High-Purity Germanium (HPGe) gamma spectrometry was used to measure the 226Ra (in secular equilibrium with 238U), 232Th and 40K activity concentration [11]. Statistical and radiological maps were produced in order to detect the relationship among values and to identify the most hazardous zone.

Moreover, since sandy natural materials were widely used in the field of civil construction, a series of radiological indices were calculated to evaluate the potential hazard connected with the building use of these materials.

### **2. Geological Setting**

The investigated river sediments come from the Ionian and Tyrrhenian seaboard of the Calabrian coastline (Figure 1), and they derive from the dismantling of the Serre Massifs and Aspromonte due to erosion processes. These Massifs belong to the Calabria-Peloritani Orogen (CPO), which is the peri-Mediterranean orogenic Alpine nappe system that comprises the whole Calabria and the north-eastern sector of Sicily [12]. It includes the Sila and Catena Costiera Massifs in northern Calabria, the Serre and Aspromonte Massifs in central and southern Calabria and the Peloritani Mountains in Sicily [12].

The Serre Massif (Central Calabria) is characterized by outcrops of granitoids and Variscan metamorphic rocks. A thin early Mesozoic sedimentary cover is present only along the Ionian margin. According to recent studies [12–14], the Serre Massif represents a nearly complete continental crust section, which is divided into lower, middle and upper crustal portions. The lower one is composed of granulitic and migmatitic complexes, whereas the middle one is made up by late-Variscan acidic granitoid rocks, such as tonalites, granodiorites and granites. Conversely, the upper crustal portion consists of high- to lowgrade metapelites (Mammola and Stilo-Pazzano complexes) [13,14].

The Aspromonte Massif (Calabria), together with the Peloritani Mountains (Sicily), constitutes the southern sector of the CPO (Calabria Peloritani Orogen). This Massif was generated by the juxtaposition of three different tectonic units. The Madonna di Polsi Unit is the lowermost unit and it is composed of greenschist to amphibolite facies metamorphic rocks showing only Alpine metamorphism. The higher-most unit is represented by the Stilo Unit, made up of phyllites, schists and paragneisses derived from the variscan metamorphism of Paleozoic sediments. The Aspromonte Unit, which is sandwiched between the Stilo and Madonna di Polsi Unit, exhibits high-grade metamorphites and late-variscan granitoids [12–14].

Paleozoic and early Mesozoic rocks of both Serre and Aspromonte Massif are covered by sedimentary deposits ranging in age from the Late Cretaceous to the Pleistocene.

### **3. Materials and Methods**

### *3.1. Samples Collection and Preparation*

Ninety samples of river sediments, around 1 kg each, were collected in nine selected sites of the Ionian and Tyrrhenian coast of Calabria for subsequent laboratory gamma spectrometric measurements. Locations, coordinates and labels of the collected samples are indicated in Figure 1 and in Table 1.

**Table 1.** Labels and Global Positioning System (GPS) locations of the sampling sites.

Samples were collected with a metal sampler and stored in labeled plastic boxes. Adequate care was taken to prevent their contamination, particularly during transportation to the laboratory. Samples were oven-dried for 24 h at 110 ◦C until moisture was completely removed and constant mass was attained. They were grounded and further homogenized by passing them through a 2 mm sieve using a sieve shaker. After, they were successfully

inserted in Marinelli hermetically sealed containers of 1 L capacity. After 40 days, the secular radioactive equilibrium between 226Ra and their daughter products was attained and samples were ready for gamma spectrometry counting.

### *3.2. Gamma Spectrometry Measurements*

Specific activity measurements of the investigated samples were performed by using a high-purity germanium (HPGe) detector (Ortec, Oak Ridge, TN, USA) for gamma spectrometry analysis.

In order to reduce statistical uncertainty, samples were counted for 70,000 seconds and spectra were analyzed in order to obtain the activity concentration of 226Ra, 232Th and 40K. The 226Ra-specific activity was calculated through a weighted average of 214Pb and 214Bi (in secular equilibrium) specific activities. The 232Th activity concentration was determined by using the 911.21 and 968.97 keV 228Ac γ-ray lines; for 40K, the evaluation was performed from its γ-line at 1460.8 keV. The experimental setup was composed by a positive biased Ortec HPGe detector (GEM) (Full Width at Half Maximum, FWHM) of 1:85 keV, peak to Compton ratio of 64:1, relative efficiency of 40% at the 1:33 MeV 60Coline), placed inside lead wells to shield the background radiation environment. Efficiency and energy calibrations were performed using a multipeak Marinelli geometry gamma source (AK-5901) of 1 L capacity, covering the energy range, 59.54–1836 keV, customized to reproduce the exact geometries of samples in a water-equivalent epoxy resin matrix [15].

The Gamma Vision (Ortec) software was used for data acquisition and analysis [16].

The activity concentration (Bq kg−1) of the investigated radionuclides was calculated using the following formula [17]:

$$\mathcal{C} = \frac{N\_E}{\varepsilon\_E t \gamma\_d \mathcal{M}} \tag{1}$$

where *NE* indicates the net area of a peak at energy *E*, *εE* and *γd* are the efficiency and yield of the photopeak at energy *E* respectively, *M* is the mass of the sample (kg) and *t* is the live time (s).

The measurement uncertainty is a combined standard one at coverage factor k = 2, taking into account the following components: uncertainty counting statistics, uncertainty in nuclear data library, uncertainty due to calibration efficiency and uncertainty about the quantity of the sample [18].

The quality of the gamma spectrometry experimental results was certified by the Italian Accreditation Body (ACCREDIA) [19]. This implies the continued verification (with annual periodicity) of the maintenance of the gamma spectrometry method performance characteristics.

In particular, with regard to the blanks, they are generally determined by acquiring an empty lead well spectrum or an ultrapure water sample in the sample holder of the desired geometry (for determination of natural-emitting gamma radionuclides) or, in the case of samples requiring a material of support (e.g., atmospheric particulate aspirated on filter), by acquiring the support of clean sampling, under the same geometric conditions as the sample. The measurement of the blanks is performed on a quarterly basis.

Regarding the precision and accuracy assessment, we proceed as follows:

For the precision, with the double test method, the repeatability of the method is verified over time. A certified reference material (also containing the radionuclides K-40, Ra-226 and Th-232) is analyzed twice, with the activity concentration of the radionuclide of interest defined as X1 (first measurement) and X2 (second measurement). The probability level p = 0.95 is considered. The following formula is applied:

$$|x\_1 - x\_2| \le \sqrt{2} \cdot s\_r \cdot t \tag{2}$$

where *t* is the Student variable and *sr* is the standard deviation of repeatability obtained in the validation phase.

For the accuracy, it is evaluated by comparing the reference (certified sample) and the measured values of the radionuclide of interest, taking into account the uncertainty. This is done using the u-test, with the following acceptability criterion:

$$u\_{test} = \frac{\text{measured} - \text{reference}}{\sqrt{u\_{\text{mcas}}^2 - u\_{ref}^2}} \le 2 \tag{3}$$

where *umeas*2 and *uref*2 are the uncertainties of the measured and reference values, respectively.

### *3.3. Radiological Indices*

3.3.1. Absorbed Gamma Dose Rate

The absorbed gamma dose rate calculation is the first major step to evaluate the health risk [20]. It was estimated as follows [21]:

$$\text{D (nGy} \cdot \text{h}^{-1}) = 0.462 \cdot \text{C}\_{\text{Ra}} + 0.604 \cdot \text{C}\_{\text{Th}} + 0.0417 \cdot \text{C}\_{\text{K}} \tag{4}$$

where CRa, CTh and CK are the mean activity concentrations (Bq kg−1) of 226Ra, 232Th and 40K in the river sediments, respectively.

### 3.3.2. The Annual Effective Dose Equivalent

Annual estimated average effective dose equivalent (AEDE) received by an individual was calculated using the following equations with an outdoor occupancy of 20% and 80% for indoors [22]:

$$\text{AEDE (outdoor) } (\text{mSv } \text{y}^{-1}) = \text{absorbed dose (nGy h}^{-1}) \times 8760 \cdot \text{h} \times 0.7 \cdot \text{Sv} \cdot \text{Gy}^{-1} \times 0.2 \times 10^{-6} \tag{5}$$

AEDE (indoor) (mSv·y−1) = absorbed dose (nGy·h−1) × 8760·h × 0.7·Sv·Gy−1 × 0.8 × 10−6 (6)

3.3.3. Radium Equivalent Activity

The activity concentrations of the radionuclides measured in the analyzed samples do not give a direct measure of the γ-radiation dose incurred by the population. To estimate the gamma radiation dose due to activity concentrations of 226Ra, 232Th and 40K in the river samples, especially when used as components of building construction, the radium equivalent activity is regularly utilized.

It is an index which describes the activities of 226Ra, 232Th and 40Kin a single activity term. The gamma energy doses from these radionuclides in river sediments are different even if they are present in the same amount. Radium equivalent activity (Raeq) estimated in Bq kg−1 is evaluated with conditions that 1 Bq kg−1 of 226Ra or 1.43 Bq kg−1 of 232Th or 0.077 Bq kg−1 of 40K produce equal gamma dose rate.

The radium equivalent activity was calculated as [23]:

$$\text{Ra}\_{\text{eq}} \text{ (Bq} \cdot \text{kg}^{-1}) = \text{C}\_{\text{Ra}} + 1.43 \text{ C}\_{\text{Th}} + 0.077 \cdot \text{C}\_{\text{K}} \tag{7}$$

### 3.3.4. Hazard Indices

External and internal radiation hazard indices were defined to limit the radiation dose to 1 mSv y−1 [24].

The external hazard index (Hex) was calculated using the given equation:

$$\rm H\_{ex} = (C\_{Ru}/370 + C\_{Th}/259 + C\_{K}/4810) \le 1\tag{8}$$

It must be lower than the unity for the radiation hazard to be negligible.

The internal hazard index (Hin) gives the internal exposure to carcinogenic radon and its short-lived progeny and it is given by the following formula:

$$\rm{H}\_{\rm in} = (\rm{C}\_{\rm Ra}/185 + \rm{C}\_{\rm Th}/259 + \rm{C}\_{\rm K}/4810) \le 1\tag{9}$$

### 3.3.5. Activity Concentration Index

Another radiation hazard, called the activity concentration index (I), has been defined by the European Commission [25] and it is given below:

$$\text{I} = \left( \text{C}\_{\text{Ra}} / \text{300} + \text{C}\_{\text{Th}} / \text{200} + \text{C}\_{\text{K}} / \text{3000} \right) \tag{10}$$

The activity concentration index is correlated with the annual dose rate due to the excess external gamma radiation caused by superficial material. Values of I ≤ 1 correspond to a criterion of 1 mSv y−1 [25]. Thus, this index should be used only as a screening tool to identify materials that might be of concern to be used as construction materials [26].

### 3.3.6. Annual Gonadal Equivalent Dose (AGED)

Some organs exhibit high sensitivity to radiation as they are highly susceptible to harmful effects of radiation. The gonads, the active bone marrow and the bone surface cells are considered to be the organs of importance due to their high radiosensitivity. The annual gonadal equivalent dose (AGED) due to the specific activities of 226Ra, 232Th and 40K is evaluated using the following formula [27]:

$$\text{AGED} \left(\mu \text{Sv} \cdot \text{y}^{-1}\right) = 3.09 \cdot \text{C}\_{\text{Ra}} + 4.18 \cdot \text{C}\_{\text{Th}} + 0.314 \cdot \text{C}\_{\text{K}} \tag{11}$$

where CRa, CTh and CK are the mean activity concentrations of 226Ra, 232Th and 40K, respectively.

### *3.4. Delineation of Radiological Map*

Surfer 10 software was employed to generate a radiological map for the study area [28]. The software uses the kriging interpolation method, used in geostatistics to predict the value of a random variable over a spatial region. Using the values of radiological measurements at all the sampling sites with different geographic coordinates, kriging allows the prediction of the radiological values throughout the regions where there are no experimental observations in the entire study area.

This method is widely employed in hydrogeology and environmental science, among others.

### *3.5. Statistical Analysis*

The XLSTAT statistical software for Windows was used for all statistical calculations [29].

With the aim to individuate the presence of the relationships among the original variables (Pearson correlation analysis), an exploratory method (Principal Component Analysis, PCA) was performed. The PCA elaboration ensures the reduction of the data dimensionality, whereas the combinations of variables identified by the PCs provide the greatest contribution to sample variability. However, before the elaboration, the logarithmic transformation of the dataset was conducted and the appropriateness of the dataset was preliminarily checked by the Kaiser–Meyer–Olkin (KMO) test and the Bartlett test (sphericity test), as already reported in our previous work [30].

### **4. Results and Discussion**

### *4.1. 226. Ra, 232Th and 40K Activity Concentration in the River Sediments, Spatial Distribution and Classification of Hazard Regions*

The activity concentrations of the three naturally occurring 226Ra, 232Th and 40K radionuclides in river sediments across the study area were determined for the sampling points reported in Table 1.

The obtained results are reported in Table 2. As shown in the Table, the highest average concentration, (31.3 ± 8.3) Bq kg−1, of 226Ra was obtained for samples from the site ID8 (Siderno–Novito river), whereas the lowest mean value, (14.5 ± 4.9) Bq kg−1, was detected for sediments from the site ID1 (Gioia Tauro–Budello river). Concerning 232Th, samples from the site ID8 are characterized by the highest average concentration, (43.5 ± 4.9) Bq kg−1, while those from the site ID2 (Gioia Tauro–Petrace river) have the lowest mean value, (20.3 ± 2.9) Bq kg−1. Sediments from the site ID1 show the highest average concentration of 40K, with (1088 ± 98) Bq kg−1, while the lowest mean concentration of (643 ± 53) Bq kg−1 was obtained for samples from the site ID 6 (Africo– Laverde river).

**Table 2.** Activity concentrations of 226Ra, 232Th and 40K in the investigated sampling sites.

Across all sampling locations, the activity concentration values of 226Ra ranged between 10.2 and 33.9 Bq kg−1, with overall mean value of 21.3 ± 6.3 Bq kg−1, that of 232Th varied from 16.5 to 48.6 Bq kg−1 and the overall average value is (30.3 ± 4.5) Bq kg−1. The 40K ranged between 528 and 1450 Bq kg−1, with an overall average value of (849 ± 79) Bq kg−1.

The average values of 226Ra, 232Th and 40K obtained in this study, for each sampling site, were compared with the world average value reported in Reference [1]. All the average values of 226Ra are lower than that of the world average (35 Bq kg−1). For 232Th, the mean values are lower than that of the world average (30 Bq kg−1) except for samples from the sites ID3 (39.5 ± 4.8 Bq kg−1), ID6 (32.5 ± 9.2 Bq kg−1), ID7 (34.5 ± 5.6 Bq kg−1) and ID8 (43.5 ± 4.9 Bq kg−1). The average concentration of 40K is higher than the 420 Bq/kg average concentration reported in Reference [1] for all the analyzed samples.

Contour maps, showing the spatial distributions of 226Ra, 232Th and 40K, are presented in Figure 2.

**Figure 2.** Contour maps showing the spatial distributions of (**a**) 226Ra, (**b**) 232Th and (**c**) 40K in the investigated region.

As can be seen in the 226Ra map (Figure 2a), there are four distinct regions. The activity concentration of 226Ra in these regions was classified as low (CRa < 15 Bq kg−1), moderate (15 Bq kg−1 < CRa < 20 Bq kg−1), high (20 Bq kg−1 < CRa < 25 Bq kg−1) and very high (CRa > 25 Bq kg−1), which are deep blue, light blue, green and yellow (and red) areas on the map. In the case of 232Th, two regions of very low and moderately low activity concentrations were identified. As seen from Figure 2b, the southern part of the study area as well as the mid-west and the northern part exhibit low activity concentrations of 232Th with typical activities concentrations less than 33 Bq kg−1. Otherwise, the middle-east zone of the study area exhibits high activity, higher than 35 Bq kg−1 . The delineated map of 40K (Figure 2c) shows two distinct zones of concentration. The blue zone in the southern part belongs to the low-concentration region (lower than 700 Bq kg−1). Whereas, other areas are classified as high-concentration zones (higher than 700 Bq kg−1).

### *4.2. Dose Assessment and Hazard Indices*

Table 3 reports the calculated average values of adsorbed and effective doses together with the hazard indices. The absorbed gamma dose rate was estimated using Equation (4). Measured values range from 52 to 79 nGy h−1, with an average value of 63 nGy h−1. This average value is near the world mean value of 60 nGy/h [1], although 66% of the analyzed samples showed values that are in excess than the world mean. The variability of the adsorbed and effective doses is attributable to the different lithologic components of the considered areas.

**Table 3.** Radiological indices in the investigated sampling sites.

Equation (5) was used to evaluate the annual outdoor effective dose due to the activities of 226Ra, 232Th and 40K in the analyzed samples. The obtained values range from 0.064 to 0.097 mSv y−1, with an average value of 0.078 mSv y−1. In a similar way, the annual indoor effective dose (Equation (6)) ranges from 0.070 to 0.323 mSv y−1, with a mean value of 0.111 mSv y−1. The total average annual effective dose is therefore of 0.189 mSv y−1. This value is lower than the worldwide average of annual effective dose (0.48 mSv y−1 [1]). The mean value is also lower than 1 mSv y−1, which is set as the maximum limit by Reference [31].

The radium equivalent activity (Raeq) was calculated using Equation (7) in order to determine the suitability of the samples when they are used as components in building materials. The obtained Raeq values vary from 52 to 79 Bq kg−1 with an average value of 63 Bq kg−1, which is lower than 370 Bq kg−1 set as the upper limit for building materials [29]. This indicates that the investigated samples may not be hazardous if used in the field of civil construction. The minimum value of radium equivalent activity was obtained

for the river sediments come from the site ID4 (Reggio Cal.–Calopinace river), while the maximum one characterizes the samples collected at the site ID8 (Siderno–Novito river).

The external and internal hazard indices, given by Equations (8) and (9), are lower than the unity for all investigated samples, thus the radiological risks can be considered negligible.

The activity concentration index, calculated with Equation (10), is again lower than the unity for all analyzed samples. Thus, these river sediments are not be of concern to be used as construction materials.

The obtained values of annual gonadal equivalent dose (AGED), calculated by Equation (11), range from 378 (site ID4) to 566 μSv y−1 (site ID8), with an average of 458 μSv y−1, for the investigated samples. These values are higher than the world average value for AGED, for sediments, which in the literature was found to be 300 μSv y−1 [32].

### *4.3. Statistical Features*

The radiological parameters obtained from the measured activity concentrations of 40K, 226Ra and 232Th were subjected to Pearson's correlation analysis, with the aim to deduce the interdependency and identify any existing relationships between the radiological indices and the concentrations of the radionuclides. This will help in making valid assessments concerning the nature and distribution of these radionuclides in the study area. The results of the Pearson's correlation coefficients are presented in Table 4. Among the negative correlations, the most interesting are those between 40K and 226Ra (−0.287) and between 40K and 232Th (−0.305), while a perfect positive correlation is found between absorbed dose rate and effective dose outdoors (1.00), and between Raeq and absorbed dose rate (1.00) and effective dose outdoors (1.00). The correlations of the three radionuclides are positive and in agreement with all the radiological parameters considered in this study at *p* < 0.01, except for correlations between 226Ra and Hin (0.739) and between 232Th and Hin (0.719), obtained at *p* < 0.05.

**Table 4.** Pearson correlation matrix among the considered variables.

Bold values refer to *p* < 0.05.

Moreover, 226Ra and 232Th are negatively correlated with the annual effective dose equivalent indoors (−0.378 and −0.281, respectively). The radiological parameters are positively and significantly correlated with each other, showing direct relationships.

In order to better define the relationships among the specific activities of 226Ra, 232Th and 40K, the PCA multivariate statistical method was carried out. The geological feature of sediments is the most important variable, neglecting the others (erosive power of rivers, solid transport, selective erosion, weathering, etc.), which is useful to evaluate the degree of relationship between the isotopic specific activities and the samples' geolocation. Results of the PCA allow to group the sampling sites in four clusters that seem to be comparable to each other (Figure 3) in terms of activity concentration of 226Ra, 232Th and 40K. Specifically, sites ID1 and ID2 show positive correlations with 40K-specific activity but negative linear relationships with 226Ra and 232Th ones, whereas sites IDs 3, 4 and 5 are negatively correlated to the activity concentrations of 226Ra, 232Th and 40K. The sites ID6 and ID7 show instead a positive correlation with values of 226Ra- and 232Th-specific activities, and exhibit a negative correlation with 40K. Finally, sites ID8 and ID9 show positive correlations of the three investigated radionuclides' activity concentrations. We believe that the differences of 226Ra-, 232Th- and 40K-specific activities among the identified clusters mainly depend on the compositional/mineralogical features of rocks from which the analyzed river sediments derive. The sediments collected along the rivers located in the Tyrrhenian side (ID1 and ID2) record the chemical fingerprint of the granitoid rocks from the Serre Massif. On the contrary, the Ionian river sediments (ID6, ID7, ID8 and ID9) reflect the chemical and mineralogical composition of the outcropping Late Cretaceous and Cenozoic terrigeous sediments. Moreover, the ID6 and ID7 rivers are characterized by a component provided by the erosion of the Aspromonte Unit, whereas clasts from the upper crustal complexes (Stilo-Pazzano and Mammola) and from granitoids supply the ID8 and ID9 rivers. Regarding the rivers on the western and southern margin of the Aspromonte Massif (ID3, ID4 and ID5), the main lithologies occurring within the drainage basins come from the Aspromonte and the Stilo Units.

**Figure 3.** Two-dimensional (2D) plots of the first two Principal Components (PCs) obtained through Principal Component Analysis (PCA) elaboration starting from 40K, 232Th and 226Ra activity concentrations (**a**). Mean concentration threedimensional (3D) map of 40K, 232Th and 226Ra into the investigated sites (**b**).

### **5. Conclusions**

The activity concentration and the distribution of natural terrestrial radionuclides 226Ra, 232Th and 40K for river sediments from selected locations of the Calabrian coastline (Southern Italy) were measured using HPGe gamma spectrometry. The obtained average values of the absorbed dose rate in air, AEDE (outdoor and indoor), Raeq, Hex and Hin, I and AGDE are comparable with average worldwide ranges and they fell within the maximum recommended values.

The background radiation maps for the study area evidence that the Tyrrhenian side shows, in general, the lowest activity concentrations of 226Ra and 232Th, whereas the Ionian zone is characterized by the highest ones. Concerning the 40K, an increasing of the specific activity was defined from south to north. Multivariate statistical analysis allowed to assess the nature and distribution of natural radionuclides in the investigated area and highlighted the relationships between chemical/mineralogical features of the local outcropping lithologies.

Data provided in this study will serve to contribute as a reference database for future researches in the area, and also as a pilot scheme to provide a comprehensive background radiation map for Italy, which is not yet available.

**Author Contributions:** Conceptualization, F.C.; methodology, F.C. and A.F.M.; validation, F.C. and G.S.; formal analysis, F.C.; investigation, F.C., A.F.M., G.S., M.D.B., M.R.F., D.R. and F.I.; resources, G.S. and G.B.; data curation, F.C. and A.F.M.; writing—original draft preparation, F.C.; supervision, G.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Applied Sciences* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/applsci

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8501-7


(5) Reform Oversight and Champions

MONRE is mandated to oversee overall water resource management but faces challenges in fulfilling its role. This is because water policy development and implementation are fragmented across national ministries with limited coordination between MONRE, other national ministries, and provincial governments [126]. To illustrate, commissions for environmental protection were established in three major river basins (Cau, Nhue-Day and Dong Nai-Sai Gon) in 2007–2009, but they have struggled to achieve environmental targets [118]. Another challenge is cooperation between national and local levels. Dual supervisory roles of central agencies and provincial people's committees, in addition to unclear reporting mechanisms, have also limited the efficiency of water governance [119,126].

Reform oversight is challenged by inadequate transboundary cooperation mechanisms. The Mekong River Commission was formed in 1995 by an agreement between Laos, Thailand, Cambodia and Vietnam, but two upstream countries of China and Myanmar are not members. As a result, decisions on building upstream hydropower plants have not fully taken into account downstream costs in relation to livelihoods and the environment [135].


**Chapter 12**

[35] Sussman DM, Schweizer KS. Microscopic theory of quiescent and deformed topologically entangled rod solutions: General formulation and relaxation after nonlinear step strain. Macromolecules. 2012;**45**:3270-3284

*Parasitology and Microbiology Research*

[36] Kalcioglu ZI, Mahmoodian R, Hu Y, Suo Z, Van Vliet KJ. From macro- to microscale poroelastic characterization of polymeric hydrogels via indentation.

Soft Matter. 2012;**8**:3393

**212**


**4. Hepatocarcinoma angiogenesis induced by DRR**

 Risk factors induced DNA damages and dysregulated DRRs are regarded as molecular events [4]. In human, risk factors for hepatocarcinoma can manifest acute and chronic DNA damage. Acute and noticeable DNA damages often lead to severe chromosome aberration and even cell death, whereas chronic DNA damages are the earliest molecular change in hepatocytes and ultimately result in hepatocarcinoma [40]. In the past decades, angiogenesis induced by dysregulation of DRR pathways may act as a vital role in the process of hepatocarcinoma. Evidence from epidemiological and clinicopathological studies has shown that higher potential of angiogenesis is in the liver of patients with chronic DNA damage and low DRR capacity [40, 102–105]. For example, Pastukh et al*.* [102] investigated the association between recruitment of DNA repair enzymes involving in BER pathway and VEGF expression via a chromatin immunoprecipitation technique. They found that hypoxia-induced reactive oxygen species (ROS) stress caused promoter base modifications targeted to hypoxic response elements (HREs) and increased VEGF expression. During this modification, 8-oxoguanine (8-oxodG, an oxidative DNA damage product) in VEGF promoter was temporally correlated with binding of human 8-oxodG glycosylase 1 (hOGG1, a BER repair enzyme), HIF-1α, redox effector factor-1, endonuclease one, and breaks in DNA strands. If 8-oxodG was decreased in the promoter region of VEGF, VEGF expression would downregulate [102]. Recent molecular epidemiological studies have further proved that genetic variants in hOGG1 genes increase hepatocarcinoma risk and modify the prognosis of this malignancy [103–105]. Collectively, these data suggest that increasing ROS like 8-oxodG resulting from low DRR capacity may promote angiogenesis.

 Studies from high HBV and HCV infection and high AFB1 exposure area also display that the degrees of DNA damages are positively associated with MVD in tumor tissues from hepatocarcinoma [20, 55, 75, 79, 82]. For example, Lu et al*.* [20] investigated the effects of XRCC4 expression in tumor tissues on clinicopathological features and prognosis of hepatocarcinoma and found that decreasing XRCC4 expression was related to low DRR capacity, causing the formation of DNA adducts and TP53M. The dysregulation of XRCC4 may promote tumor proliferation and increase MVD. Several other studies further show that the low DRR capacity resulting from significant mutations in coding region of DNA repair genes (such as XRCC4, XRCC1, XPC, XPD, and XRCC7) increases MVD (**Table 4**) [21, 40, 52, 55, 59, 61, 62, 79, 80, 82]. Results from Lu et al*.* [20] and our studies [61, 62] showed that genetic alterations in the coding regions of XRCC4 gene (including Ala to Ser at codon 247 and Thr to Ala at codon 56) can decrease levels of XRCC4 protein expression and cause increasing amount of AFB1-DNA adducts and mutative frequency of TP53 gene in tissues with hepatocarcinoma. They also found that the amount of AFB1-induced DNA adducts, including 8,9-dihydro-8-(N7 -guanyl)-9 hydroxy-AFB1 (AFB1-N7 -Gua) and formamidopyridine AFB1 adduct (AFB1-FAPy), was positively associated with the number of microvessels (a biomarker for angiogenesis). Results from our studies [79, 106, 107] furthermore displayed that three low DNA repair markers related to AFB1, including tumor risk, TP53M frequency, and AFB1-FAPy adduct amount, were significantly correlated with the number of microvessels in liver tissues. These individuals with high AFB1-FAPy adduct level in liver tissues had an increasing risk of high MVD than those low adduct level (OR = 1.68, 95% CI = 1.45–2.87) [106]. Liu et al*.* [108] and Wang et al*.* [109] further proved that the upregulation of microRNA-429 and microRNA-24 expression in tissues with hepatocarcinoma not only increased the amount of AFB1-DNA adducts

 and the number of microvessels but also grew tumor metastasis risk via vessels and shorted patients' survival. Recent evidence has shown that microRNA-24/ microRNA-429 can modify the capacity of DDR via controlling Nbs1 (a regulator of DRR) [110, 111] and angiogenesis via regulating the crosstalk between the pro-contractile transforming growth factor-β/bone morphogenetic protein (TGF-β/ BMP) signal (inducing a quiescent 'contractile' phenotype) and the pro-synthetic platelet-derived growth factor (PDGF) signal (causing a proliferative 'synthetic' phenotype) [112, 113]. This suggests that microRNA-24/microRNA-429 may play an important regulative role between DRR capacity and angiogenesis. Taken together, this evidence proves that low DRR-induced MVD augmentation is regulated by the amount of DNA damage.

 Evidence from in vitro and in vivo studies further shows that dysregulation of DRRs and signaling to cell cycle checkpoints (CCCs) may modify hepatocarcinoma angiogenesis. CCCs involving in DRRs mainly encompass G1/S and G2/M checkpoint [114]. During G1/S checkpoint, both ATR and ATM act as central activators for DRR via inducing the phosphorylation of p53 protein which can activate p21 (a Cdk inhibitor). ATM/TP53/P21 pathway also plays an important function controlling G2/M procession [114]. The dysregulation of these factors and signal pathways can change the status of angiogenesis [115–119]. For example, Qin et al*.* [115] found that E2F1, an important cell cycle regulator, can modify angiogenesis via controlling VEGF expression by p53-dependent way. In this control model, deficient phenotype of E2F1 will result in VEGF overexpression, while its positive phenotype decreases VEGF expression [115]. Factors controlling cell shape and cytosol can regulate the cycle of vessel endothelial cells and angiogenesis [116, 117]. In mice model with the deficiency of BCL-2 (an important regulatory factor in DDRs), cells featured increasing DNA damage [118]; the inhibition of BCL-2 will result in the arrest of cells in S phrase and suppression of tumor angiogenesis [119]. In an integrated genomic study (including 5 hepatocarcinoma patients with hepatitis D visus [HDV] and 7 HDV-positive cirrhosis cases), Diaz et al*.* [120] investigated the association between HDV-related hepatocarcinoma and potential signal pathways involved in DNA damage and repair and cell cycle and found significant interactions of DDR/cell cycle-related genes, such as BRCA1, BARD1, CDK1, CDKN2C, CCNA2, CCNB1, CCNE2, GSK3B, H2AFX, MSH2, NPM1, PRKDC, and TOP2A. Results from the t-SNP (*t*-distributed stochastic neighbor embedding analyses) further exhibited that HUS1, BRCA1, BARD1, GADD45, DNA-damage-induced 14-3-3σ, and MSH2 gene involving in DRRs valuably scored with regulatory genes (such as ATM, TP53, NO, and epidermal growth factor), which involve in G2/M checkpoint and angiogenesis [120]. The dysregulation of HUS1 and corresponding genotoxin-activated checkpoint complex (also termed as Rad9- Rad1-Hus1complex) will cause abnormal DRR capacity and cell cycle in response to DNA damage and promote the alteration of hematogenous metastatic phenotype for hepatocarcinoma [121, 122]. The genetic alterations and abnormal expression of BRCA1 and GADD45 (two important regulatory factors in DRR and apoptosis pathways) in hepatocytes can also change TP53-dependent CCCs and VEGF expression [123, 124]. Altogether, these studies have proved that the dysregulation of DDRs can cause the abnormal regulation of CCCs and change the status of hepatocarcinoma angiogenesis.

Detailed molecular mechanisms of DRR dysregulation promoting hepatocarcinoma angiogenesis have still not been fully understood. Several possible pathways may play some important roles. First, DNA damage agents induce NO synthase and increase the expression of VEGF and HGF [125, 126]. Second, DNA damage agents like AFB1 cause the mutations of such genes as TP53, ras, and DNA repair genes. Activation of oncogenes and inactivation of tumor suppression genes and

 DNA repair genes lead to uncontrolled expression of genes involving in angiogenesis such as VEGF and Ang-1/2 [5, 6]. Third, genetic alterations in DRR pathways may alter the microenvironment of tumor and promote angiogenesis [127–129]. Fourth, the abnormal DRRs may accelerate the accumulation of DNA damages and trigger the dysregulation of angiogenesis-related genes and the progression of hepatocarcinoma. Finally, some metabolic products (such as AFBO) or nucleotide sequences (HBx) of DNA agents can bind to genomic DNA of hepatocytes and may increase the activation of VEGF HREs [22, 40, 41, 45]. Taken together, under the conditions of low DRR capacity and/or chronic risk factors, DNA damages will accumulate in hepatocytes and ultimately induce hepatocarcinogenesis and tumor angiogenesis.


**1. Introduction**

In recent decades, environmental policies have been oriented towards reducing energetic dependence on fossil fuels, leading to renewed and increased investments in the hydroelectric sector; such investments have been primarily directed to optimizing both new designs and existing installations, thus allowing for more efficient plants and reduced payback times. Kaplan hydro turbines have been adopted for a long time to deal with efficient energy production in the range of high specific speeds [1]. The use of movable vaned distributors and runner blades, in fact, allows to limit incidence flow losses and achieve high conversion efficiency over a wide range of flow rates, from 40% to 120% of the nominal value. The main drawback of Kaplan turbines is that their size increases as the flow rate gets higher, and the cost attributed to civil works becomes more and more substantial.

In these machines, the draft tube is an essential component as it makes it possible to recover a significant fraction of the kinetic energy leaving the runner by conversion into static pressure. This, in turn, leads to creating an additional suction head downstream of the runner, which allows to increase the effective head that the runner is able to deliver [2]. In addition, the draft tube plays a fundamental role in determining turbine efficiency since the height at which the tube is installed contributes to a large percentage of the total net head that can be recovered [3].

Moreover, it is one of the most challenging parts to analyze from a fluid flow perspective due to the interaction of many complex flow features, such as unsteadiness, turbulence, separation, streamline curvature, secondary flow, swirl, and vortex breakdown [4].

Studies involving Kaplan turbine components made use of both steady and unsteady Reynolds-Averaged Navier-Stokes (RANS), the Scale-Adaptive Simulation (SAS) and hybrid RANS-LES (Large Eddy Simulations) models, such as the Zonal Large Eddy Simulation (ZLES) or the Detached Eddy Simulations (DES). Jošt reported an improvement in characteristic parameter predictions for a Kaplan turbine using unsteady models and ZLES [5]. Marjavaara [4] performed steady and unsteady CFD simulations of a draft tube using RANS and DES. The characteristic quantities and the local flow features obtained were highly dependent on the model and the inflow boundary conditions. These models can accurately predict turbulent phenomena in the draft tube but require very large computational resources. As such, their adoption in optimization studies seems, today, still infeasible, while validated RANS approaches offer a very good compromise between accuracy and computational effort [3,5]. Li [6] used steady RANS to study the hydraulic performance of a Kaplan turbine's components. Brijkishore [7] evaluated the influence of runner solidity using the k-Ω SST turbulence model in steady computations. Liu [8] carried out transient analyses on the runner lifting-up during load rejection using the same RANS model.

Past studies on draft tube optimizations were limited in terms of representativeness. In fact, while the draft tube can be acceptably optimized as an isolated component around its design point, therefore neglecting the interference effects with the turbine runner, the performances of optimal configurations have not been post-checked using the entire turbine model for a posteriori validation of the complete turbine installation. Lyutov [9] described a coupled optimization of runner and draft tubes, where their geometries were varied simultaneously, allowing to achieve a 0.3% higher efficiency gain compared to single-runner optimization. Ciocan [10] presented a draft tube optimization from part to full load, using a parameterized velocity field at the tube entry based on a swirl-free velocity profile for the runner. Eisinger [11] reported an automatic draft tube optimization, but no coupling of the inflow conditions with the attached runner was present. Puente [12] employed RANS simulations with the k-ε turbulence model for the same purpose, using normalized velocity profiles at the duct inlet, derived from previous turbine simulations. Sosa [13] analyzed the influence of draft tube modification using a steady computational model with a single-channel guide vane and runner connected to the tube via a mixing plane. The geometric modifications based on a two-parameter description were meant to represent erosion or material addition in the wall. Schiffer [14] performed steady-state simulations of a complete model, including distributors, runners, and a draft tube, at several operating points. The tube was then manually adapted to optimize the turbine's efficiency while keeping construction costs low.

In the present paper, a fully 3D Kaplan turbine RANS CFD model was first built and validated against available experimental data in terms of net head, flow rate, efficiency, and mechanical power.

Following this, the pressure and velocity distribution at the draft tube inlet were extracted from previous computations. Finally, a CFD model that contained only the draft tube and its outlet extension was created and used in the optimization process. Unlike previous studies, the optimization was carried out based on two-objectives which are of prominent interest to draft tube designers, i.e., its pressure recovery and its total pressure loss coefficient. Also, a noteworthy improvement with respect to past studies relies on the verification of the obtained results: in fact, once the optimal draft tube geometries had been found, they were tested back using the full machine CFD model to find out their impact on the turbine's efficiency and to analyze the sources of draft tube losses in detail.

#### **2. Materials and Methods**

#### *2.1. Draft Tube Performance*

The datum, or baseline, draft tube was connected to a Kaplan turbine runner installed in a 576 kW hydraulic power plant located in the Adda river (Italy), as depicted in Figure 1a. The original geometry adopted for the elbow draft tube is schematically represented in Figure 1b.

**Figure 1.** Original power plant: (**a**) overall external view; (**b**) baseline elbow-draft tube geometry.

As mentioned above, the draft tube behaves as a diffuser acting to recover part of the kinetic energy from Section 2 to 3, according to its prescribed area distribution along the equivalent centerline.

Two coefficients are used here to quantify the performance of a draft tube, namely a pressure recovery factor (*Cp*) and a drag coefficient (*Cd*), respectively defined as:

$$\mathcal{C}\_p = \frac{p\_3 - p\_2}{\frac{1}{2}\rho c\_2^2}; \; \mathcal{C}\_d = \frac{p\_{\text{tot}-2} - p\_{\text{tot}-3}}{\frac{1}{2}\rho c\_2^2} \tag{1}$$

where *p* is the static pressure and *ptot* the total pressure; subscripts 2 and 3 refer to the station downstream of the runner and on the draft tube outlet, respectively; ρ is the water density; *c* is the absolute velocity.

The values of these coefficients are strictly related to the area distribution and are notoriously conflicting in such a way that the more intense the diffusion is (high *Cp*), the bigger the head loss is (high *Cd*) that might occur as a result of flow separation, secondary losses, etc. Therefore, a good draft tube features the highest possible *Cp* along with the minimum achievable *Cd*. These two coefficients were used as objective functions in the multi-objective optimization study described later.

#### *2.2. CFD Model Setup and Validation*

The steady-state RANS CFD model of the complete turbine adopted included complete spiral casing with stay vanes, guide vanes, and runner and draft tubes (Figure 2). The final grids were chosen based on a mesh sensitivity analysis and their details and statistics are provided in Table 1. The influence of mesh size on the draft tube performance parameters is shown in Figure 3, where a limited sensitivity to the grid size can be observed, starting from the third level, having 1.8 M cells. As for boundary conditions, a total pressure was defined at the inlet of the spiral casing and a static pressure equal to 0 Pa (relative to the atmospheric pressure value of 1 atm) was set at the outlet of the draft tube.

The total pressure at the inlet boundary corresponded to the total energy that the plant could process and was calculated from the net head and the flow rate that the turbine was subjected to. The flow solver was Ansys CFX 19.2. Since the runner was not modelled as a rotating mesh, a multiple-reference-frame approach was chosen, and a mixing-plane condition was set both at the distributor–runner and runner–draft tube interfaces. Among all the turbulence models (all the variants of both k-ε and k-ω available in Ansys CFX 19.2 were actually tested), it was observed that the k-ω SST [15] provided the most accurate predictions, but it was also the one that required the highest computational effort; for this reason, the k-ω SST was used in off-loop analyses on the complete machine model, while the standard k-ε was used in the optimization due to its cheaper requirement. A physical timescale corresponding to 1◦ of runner rotation and a high-resolution scheme were set.

In order to assure improved numerical stability, 50 iterations were run at first using first-order schemes. These results were adopted as initial values for the high-resolution scheme.

**Figure 2.** Turbine computational domains: (**a**) close up of the runner domain in meridional view; (**b**) 3D view of the overall domain; (**c**) close up of the computational mesh in the region close to the runner/draft tube interface; (**d**) top view of the overall domain; (**e**,**f**) details of the draft tube mesh.

**Table 1.** Details of the grids used in the subdomains.

**Figure 3.** Mesh independency analysis for the draft tube.

During the convergence runs, several variables of interest were monitored, such as flow rate, efficiency, static pressure at draft tube inlet, and mechanical power output. Final results were available after 3000 iterations on average and almost 72 h of CPU calculations. All calculations were performed on a server equipped with Intel® Xeon® CPU x5650 processors using parallel solutions on 10 multiple cores. Convergence was assessed when all variables of interest showed a variation lower than 0.08%.

Results from both experiments and simulations are presented in Figure 4. All the experimental data (net head, flow rate, efficiency, and mechanical power) were taken in compliance with the European IEC EN 60041:1991 standard [16] and included the appropriate uncertainty bandwidth. Power and efficiency data from CFD were calculated as follows:

$$P = M \cdot \omega; \; \eta = \frac{P}{\rho Q g H} \tag{2}$$

where *M* is the total torque on the runner blades, ω is the rotational speed, *Q* is the flow rate, and *H* is the net head. As can be observed in Figure 4, the power and efficiency data were within the experimental measurement uncertainty range. The steady RANS CFD model set up appeared, therefore, to be able to reproduce the machine performance with sufficient accuracy and acceptable computational time. From this validated model, the runner–draft tube interface conditions, that is, the distributions of the *u*, *v*, *w* (respectively tangential, absolute, and relative velocity components), were extracted to be used as inlet boundary conditions during the optimization process.

**Figure 4.** Validation results: (**a**) power data; (**b**) efficiency.

#### *2.3. Draft Tube Geometry and Parameterization*

The draft tube geometry of Kapan turbines is complicated by the di fferent shapes of its inlet and outlet sections. The first is circular due to the interface with the runner, while the second is often rectangular. Moreover, the draft tube features a 90◦ elbow to minimize excavation costs and to improve powerhouse compactness. As a result, a large number of design parameters is necessary to provide successful geometry parameterization. In the present paper, the draft tube was parameterized using a mean line and a number of cross-sections stacked along it.

The mean line was designed to be composed by a first straight section, made up of two segments, one related to the divergent cone, and one to the cylindrical section; a second curvilinear line was related to the elbow shape and controlled by a three-point Bézier polygon; finally, a segment was connected to the last section of the draft tube.

The first straight part of the mean line was described by two parameters—the two segments' lengths. For the Bézier polygon, four additional variables were introduced—the coordinates of the polygon-points (the coordinates of the first Bézier polygon point were derived from the first straight part). The last parameter was the length of the last straight segment. In total, 7 variables were involved in the parameterization of the mean line. The area distribution was managed using a 4th-degree polynomial as a function of the mean line coordinate *x*:

$$A(\mathbf{x}) = P\mathbf{1} \times \left(P\mathbf{2} \times a\mathbf{x}^4 + b\mathbf{x}^3 + c\mathbf{x}^2 + d\mathbf{x}\right) + \varepsilon \tag{3}$$

In Figure 5, the area distribution of the original geometry is represented using dots along with the polynomial interpolation. From the original geometry, it was possible to obtain the polynomial coe fficients (*a*, *b*, *c*, *d*, *e*) and use P1 and P2 as decision variables.

**Figure 5.** Draft tube geometry definition: (**a**) area distribution; (**b**) shape of a generic cross section.

Overall, the draft tube was parameterized using 9 decision variables. In order to avoid the need for further parameters to account for the shape change from inlet to exit sections, the geometric similarity with the baseline tube was preserved, keeping the geometric ratios R1/L, R2/L, and H/L for each section along the mean line equal to those of the datum.

#### *2.4. Optimization: Problem Formulation and Tools*

As previously mentioned, the optimization study was conducted on the draft tube alone by taking the inlet boundary conditions from the full machine model at the runner–draft tube interface. As a result, the flow domain included the runner discharge cone along with the runner tube interface

(which form the inlet boundary), the first divergent part, the elbow, the last divergent part, and the draft tube prolongation leading to the outlet section.

The optimization was multi-objective, with the purpose of maximizing *Cp* and minimizing *Cd* at the same time using a Pareto approach. As for the constraints, the optimization had to lead to solutions that were interchangeable with the baseline case currently installed without any modifications regarding the runner. As a result, during the entire procedure, the baseline inlet and runner cone (both belonging to the runner) remained fixed. A further constraint was related to the global draft tube depth which a ffects the cost of civil works. The outlet depth could vary in all directions but with the limit of always remaining with a certain margin under the free water surface of the tailrace.

A response surface methodology for the optimization tool was used [17]. The first step was a design space exploration using design of experiments (DOE) in order to understand how the design parameters were related to each other and which were the most significant ones. A "central composite design" (CCD) DOE technique was adopted, which provided a screening set of 150 samples to determine the overall trends of the metamodel. The latter was built using the genetic aggregation method [18] and the Pareto optimal solutions were finally obtained (Figure 6a). A check-in-the-loop of optimal solutions was necessary to evaluate the response surface goodness using CFD (Figure 6b); the best-so-far optimum individuals from the response surface were taken and verified using a dedicated CFD simulation, then inserted back in the DOE, which, in turn, updated the response surface from which a new Pareto front was obtained. The loop was kept running until a stable and invariant Pareto optimal set was determined.

**Figure 6.** Optimization procedure: (**a**) general loop; (**b**) response surface update workflow.

It is worth mentioning that during the loop iterations, the response surface method allowed for an adaptive refinement in the search regions around maximum *CP* and minimum *CD* so as to make it possible to have a wider and more uniform Pareto front.

#### **3. Results and Discussion**

The final Pareto front is depicted in Figure 7, from which a subset of relevant configurations could be extracted: the baseline solution is marked with a black diamond, while the optimum *CD*, the optimum *CP*, and a good compromise between the *CD* and *CP* are marked using a red triangle, a yellow circle, and a green square, respectively. Dominated individuals are indicated as other points in the figure. The performance data of optimized individuals are given in Table 2. The solution corresponding to the minimum *CD* (Candidate 1) featured a reduction count Δ *CD* = −0.0254 (i.e., −2.54 percentage point reduction), while the one corresponding to the maximum *CP* (Candidate 3) showed Δ *CP* = +0.0574 (i.e., +5.7 percentage point improvement).

**Figure 7.** Optimization results: Pareto solutions.

#### *3.1. Optimal Draft Tube Configurations*

In Figure 8, the area distribution along the mean line coordinate for the Pareto optimal subset of solutions compared to the baseline one is shown. It can be seen that optimal solutions exhibited a much lower increase in area distribution from the beginning of the draft the tube (mean line coordinate > 0.05, see Figure 8b); this is beneficial to limit adverse phenomena during flow diffusion in this zone as confirmed, also, by the distribution of Cd versus the mean line coordinate (Figure 10a), particularly visible in the case of Candidate 1.

**Figure 8.** Optimization results: (**a**) area distribution (**b**) close up in the range *x* ∈ [0–0.3].

In Figure 9, a detailed comparison between the optimal solutions and the baseline is presented. For the first draft tube part (until the 20% of the mean line), very similar characteristics can be observed, while the most remarkable di fferences are visible in the elbow and in the last part.

**Figure 9.** Comparison among geometries in meridional view: Baseline (grey); Candidate 1 (red); Candidate 2 (green); Candidate 3 (yellow).

Candidate 1 had a very similar geometry compared to the baseline, while Candidate 2 featured an elbow with a much higher radius of curvature. Candidate 3, on the other hand, had much larger dimensions than the baseline and this has a major impact on the excavation costs. As it will be seen later, Candidate 2 was the best in performance and this is due to the shape of the elbow, a fact that emphasizes the importance of producing a good design in the region where a big part of the di ffusion is realized by the draft tube.

#### *3.2. Post-Check Validation on Full Machine and Result Transposition*

A number of Pareto optimal draft tube geometries were obtained from the optimization and were subsequently analyzed, a posteriori, by connecting them into the CFD model of the entire machine in order to assess their influence on critical plant characteristics, such as mechanical power and e fficiency. In Figure 10, both the *Cd* and *Cp* of the optimal draft tube geometries computed using the full CFD model of the turbine are plotted as functions of the mean line coordinate.

**Figure 10.** Post-check results on draft tube performance using full machine model: (**a**) *Cd* distribution; (**b**) *Cp* distribution.

The CFD analyses performed on the full machine returned di fferent values for the turbine performance figures, especially for the net head and mass flow rate (Table 3). For a better comparison, a result transposition was carried out according to the European IEC EN 60041:1991 standard [16]. In particular, the mass flow rate and the power output data were transposed using the baseline net head. This was acceptable since the condition 0.99 < √H/ √H' < 1.01 on head H' compared to baseline H was met for each candidate. Optimal candidate solutions led to greater power and e fficiency values compared to the baseline. In particular, Candidate 2 featured an increment of 1.8% on the produced power and 1% in e fficiency ratio.

**Table 3.** Full machine CFD results (transposed based on baseline CFD results).

**Figure 11.** Comparison of geometries (left column), and streamlines colored by velocity (right column): (**a**,**a'**) Baseline; (**b**,**b'**) Candidate 1; (**c**,**c'**) Candidate 2; (**d**,**d'**) Candidate 3.

From a geometrical point of view (Figure 9), Candidates 1 and 2 were not very different from the baseline, although they clearly exhibited an elbow shape having a larger radius of curvature, which eventually caused the exit stations of both the elbows and the draft tubes to be placed further downstream with respect to the baseline. Additionally, Candidate 3 featured the same global tendencies but with a much longer draft tube, which, eventually, would lead to higher installation costs. In Figure 11 a 3D view of the several solutions is presented, with a comparison of the several streamlines. It can be seen that candidate 1 has the most uniform flow, a fact confirmed by his *Cd* value (Figure 9).

From a fluid dynamic standpoint, it is apparent that the largest part of the pressure recovery took place in the first part of the diffuser (Figure 10b), which confirms the importance of having a very accurate design of the region upstream of the elbow.

As a matter of fact, approximately 70% of the pressure recovery is realized in the first 20% length of the diffuser. This also held true for the total pressure losses which were primarily generated in this region (see C*d* behavior in Figure 10, where the largest slope is evidenced in the C*d* distribution in the first 20% length of the baseline geometry); localized losses accumulated all along the draft tube mean line, finally leading to an overall large C*d* value in the baseline. Such an observation confirms what was previously pointed out in [5] and [13]. Furthermore, by looking at the total pressure distribution on consecutive station cuts along the draft tube (Figure 12, right column), it is evident that the rotational region at the lowest total pressure in the core of the discharge flow was of a much lower intensity in the optimized solutions compared to the baseline, despite the latter having a separation wall in the middle of the diffuser. In fact, the analyses carried out showed that the vortex rope greatly influences the draft tube's performance, creating instabilities and turbulent zones in the flow field.

**Figure 12.** Comparison of relative total pressure distributions: (**a**) Baseline; (**b**) Candidate 1; (**c**) Candidate 2; (**d**) Candidate 3.

## **4. Conclusions**

A validated CFD model of a draft tube of a Kaplan turbine was successfully implemented and used for a multi-objective optimization based on the construction of a response surface. The shape optimization involved only the draft tube domain, where a one-way coupling with the runner was included in the form of a prescribed velocity distribution extracted from the full-component baseline model.

Two Pareto optimal designs have been extracted from the final set which outperformed the baseline. Regarding the performance of the draft tube alone, the optimum CD featured a reduction count of –0.0254 (–2.54 percentage point reduction), while the optimum CP showed an incremental count of +0.0574 (i.e., +5.7 percentage point improvement) with respect to the baseline geometry.

A post-check carried out including the optimized draft tubes in the overall turbine confirmed the optimization trends, although the best improvements in terms of delivered hydraulic power—in the case where the overall turbine was considered—have been registered, including a third type of diffuser which exhibited an increment of 1.8% in the power produced.

The three best candidates overall were slightly deeper than the baseline in terms of excavation dimensions, so a cost/benefit assessment will be needed to determine which is the most appropriate draft tube for a given plant. In addition, future work will be related to extending the investigation to more operating points by performing a multi-objective multi-point optimization, where the Pareto solutions will account for the different conditions at the tube inflow along the turbine operating line.

**Author Contributions:** Conceptualization, E.B. and R.O.; methodology, R.O. and E.B.; software, R.O.; validation, R.O., M.M. and R.B.; formal analysis and writing—original draft preparation, R.O.; writing—review and editing, E.B.; supervision, E.B.; review and editing, A.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.