**Green Roofs as E**ff**ective Tools for Improving the Indoor Comfort Levels of Buildings—An Application to a Case Study in Sicily**

**Laura Cirrincione 1,2,\*, Maria La Gennusa <sup>1</sup> , Giorgia Peri <sup>1</sup> , Gianfranco Rizzo <sup>1</sup> , Gianluca Scaccianoce 1,3, Giancarlo Sorrentino <sup>1</sup> and Simona Aprile <sup>4</sup>**


Received: 27 December 2019; Accepted: 22 January 2020; Published: 29 January 2020

**Abstract:** In the line of pursuing better energy efficiency in human activities that would result in a more sustainable utilization of resources, the building sector plays a relevant role, being responsible for almost 40% of both energy consumption and the release of pollutant substances in the atmosphere. For this purpose, techniques aimed at improving the energy performances of buildings' envelopes are of paramount importance. Among them, green roofs are becoming increasingly popular due to their capability of reducing the (electric) energy needs for (summer) climatization of buildings, hence also positively affecting the indoor comfort levels for the occupants. Clearly, reliable tools for the modelling of these envelope components are needed, requiring the availability of suitable field data. Starting with the results of a case study designed to estimate how the adoption of green roofs on a Sicilian building could positively affect its energy performance, this paper shows the impact of this technology on indoor comfort and energy consumption, as well as on the reduction of direct and indirect CO<sup>2</sup> emissions related to the climatization of the building. Specifically, the ceiling surface temperatures of some rooms located underneath six different types of green roofs were monitored. Subsequently, the obtained data were used as input for one of the most widely used simulation models, i.e., EnergyPlus, to evaluate the indoor comfort levels and the achievable energy demand savings of the building involved. From these field analyses, green roofs were shown to contribute to the mitigation of the indoor air temperatures, thus producing an improvement of the comfort conditions, especially in summer conditions, despite some worsening during transition periods seeming to arise.

**Keywords:** innovative envelope; building components; green roofs; indoor comfort; energy consumption; building modelling; simulation models

#### **1. Introduction**

The reduction of energy consumption and the related decrease of greenhouse gases emissions represent important aspects to which much attention has been paid at global, European and countries levels, especially with regard to the building sector.

Worldwide, energy consumption in the building sector is responsible for 36% of total energy use (corresponding to a 39% of energy-related CO<sup>2</sup> emissions) [1,2], while at the European level, the energy consumption in the same sector accounts for a share of the total energy comprised between 25% and 40% (corresponding to about 35% of CO<sup>2</sup> emissions throughout Europe) [3–5].

From this perspective, various strategies have been implemented. At the global level, the UN 2030 Agenda for Sustainable Development, along with the 17 Sustainable Development Goals (SDGs) [6], need to be mentioned.

At the European scale, the EU has been very committed to this issue by setting the well known ambitious targets for 2020 ("climate and energy package") [7], and even more so for 2030 ("climate and energy framework") [8,9] and 2050 ("long-term strategy") [10,11]. Other relevant goals have been set out in the seventh Environment Action Program (EAP) [12] aimed at decarbonizing and making more sustainable European cities. Among European standards and regulations issued on this matter, the EPBD Directive and its recast must be cited [13–15].

Italy's National Energy Strategy 2017 [16] lays down the actions to be achieved by 2030, in accordance with the long-term scenario drawn up in the EU Energy Roadmap 2050, which translate to a reduction of emissions by at least 80% from their 1990 levels.

However, despite these standards and regulations being in force, in recent years, the energy consumption in the building sector has increased, particularly in Italy [17]. That is why more effort in promoting actions and finding new strategies to improve energy savings and efficiency are necessary [18].

Generally speaking, apart from all the design strategies typical of the principles of bioclimatic architecture (such as, for instance, space organisation, wall-window-ratio, orientation, thermal mass, operation management [19,20]), more relevant energy savings achievable in buildings can be attributed to two main categories of components: technical plants (HVAC system) and the building envelope, which have a synergistic relationship. In fact, a reduction in energy consumption related to the HVAC consists in the use of active systems which entail further energy consumption. As regards the building envelope, passive systems (not energy depending) can be used, which allow to actually obtain a reduction of the energy consumption (and at the same time, to also save on the use and the size of the HVAC system). Clearly, the occupants' behaviours and attitudes might also significantly influence energy saving, as demonstrated, for instance, in [21–24]. Starting from the above considerations, in this work, it was decided to pay attention to the use of a passive system to be applied to the building envelope, that is, green roofs equipped with different vegetation types.

Among the passive systems, green roofs are becoming more and more popular due to their capability of reducing the energy needs for the climatization of buildings [25–27], especially for cooling purposes [28–30].

At the same time, vegetated roofs also have a positive impact on the outdoor urban environment in terms of regenerative sustainability, allowing to induce various environmental benefits [26,31], such as reducing air pollution [32,33], mitigating noise [34,35], improving the management of runoff water [32,36,37], easing the urban heat island (UHI) effects [38–40], and increasing the urban biodiversity [41,42]. Moreover, the European Union is evaluating the possibility of including criteria specifically referring to green coverings within the EU Ecolabel scheme for buildings [43].

In addition to experimental studies [44–46], the effect on the built environment of vegetated roofs in diverse climates has also been investigated from analytical [47–49] and modelling points of view [50–52] over the years. In particular, the relevant parameters for energy modelling of green roofs have been explored in the literature, particularly referring to the role played by leaves and solar radiation in the thermal exchanges between vegetated layers and the surrounding environment [53].

The reported literature indicates that green roofs represent very promising building components, also in the Mediterranean context, as demonstrates the incremental number of studies and analyses carried out in recent years concerning both the experimental [44,54] and the simulating approach [47,55].

Other studies have underlined that additional issues would probably need more attention regarding plants growing on the roof, especially their influence on the thermal performance of green roofs [56,57], and the influence of the evapotranspiration component on the green roof heat and mass transmission [50,58].

Taking into account the studies reported above, it is evident how green roofs can have a strong impact in attenuating the average radiant temperature on building roofs [44,59,60]. This capability of acting as thermal insulation positively influences the indoor comfort conditions for the occupants of the rooms sited under the roof [27,45,61,62]: this aspect has always been critical in the design phase of a building envelope.

Kuan-Teng Lei et al. [63], by means of a field experiment performed in a school building in Taipei, have developed a finite element analysis model for the improvement of indoor thermal comfort in the presence of extensive green roofs. The researchers found a decrease of the indoor temperature up to 4 ◦C, compared to bare roofs. Costa Junior et al. [64], through an experimental analysis conducted in the city of Recife, Pernambuco (Brazil), compared the performance of four roofs made up of chanana green roof (Turnera subulata), daisy green roof (Sphagneticola trilobata), parsley green roof (Ipomoea asarifolia), and fiber cement tile. Through the comparison, the index of discomfort (ID), effective temperature (ET) and the human comfort index (HCI) were calculated. The three vegetated options mitigated both the internal air temperature with a reduction of 0.71◦C, 0.19◦C and 0.35◦C, respectively and the internal surface temperature with a reduction of 1.5◦C, 0.8◦C and 0.8◦C, respectively, compared to the fiber cement tile-made roof. Di Giuseppe and Orazio [65] experimentally analysed the effect of cool and green roofs compared to traditional ones in a Nearly Zero Energy Building, on the internal comfort and the air temperatures of the surrounding environment. The outcomes, on one hand, confirmed the effectiveness of green and cool roofs for the mitigation of the Urban Heat Island effect, and on the other hand, indicated the little effectiveness of high-albedo materials on roofing systems with a very low U-value for internal comfort.

Furthermore, the impact of green walls on thermal comfort have been compared to that of green roofs. For instance, Malys et al. [66], using the SOLENE-microclimat tool, compared the effect caused by different 'greening strategies' on buildings' energy consumption and indoor comfort in the summer season. The outcomes of the investigation indicate that, while green roofs seemingly mainly affect the upper floor, green walls directly affect the indoor comfort throughout the entire building.

To help to provide a contribution to this important and often overlooked matter, the aim of this paper was to assess the influence that green roofs have on the indoor thermal comfort levels, particularly considering the indoor radiative heat exchanges. For this purpose, a case study was conducted to estimate how the adoption of the proposed interventions could impact the indoor thermal comfort and the energy consumption of a building and contribute to the reduction of the direct and indirect CO<sup>2</sup> emissions. In particular, the ceiling temperatures of some rooms located underneath six different types of green roofs were monitored. The choice of detecting this parameter resides on the circumstance that the ceiling internal surface's temperature is a relevant component of the mean radiant temperature of the room that, in turn, greatly affects the value of the indoor parameter PMV [67]. Subsequently, the obtained data were used as input data for one of the most widely used simulation models (EnergyPlus [68]) to evaluate the indoor comfort levels and the achievable energy demand savings of the involved building.

Three scenarios were adopted. Scenarios #1 and #2 refer to a building equipped with a green roof; Scenarios #3 refers to the pre-existing roof. The simulation of Scenario #2 is made by means of the Energy Plus code, through its resident routine; on the other hand, Scenario #1 is modelled by imposing, for the indoor temperatures of the ceiling, the experimental data detected on the site. The aim of this approach was to compare the PMVs obtained from Energy Plus with those calculated on the basis of the monitored experimental data, i.e., the indoor ceiling surface temperatures.

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

The presented study is part of a joint research project between the University of Palermo and the "Consiglio per la ricerca in agricoltura e l'analisi dell'economia agraria - CREA", operating in Sicily. To accomplish the task mentioned in the Introduction, a mixed approach, partly modelling and partly experimental, was used in the work.

At the same time, the impact that green roofs have on the energy consumption of a building was evaluated. In addition, the estimation of the achievable savings in direct and indirect CO<sup>2</sup> emissions due to the use of such building component is reported as well.

#### *2.1. Description of the Experimental Site*

The installation of the experimental green roof was settled by the CREA Research Center and the University of Palermo with the support of a building materials enterprise, on the roof of a one-storey detached house (Figure 1) owned by CREA and sited in Bagheria, a Sicilian town near Palermo (Southern Italy).

**Figure 1.** Site of the installation of the experimental field (source: Google-Earth).

To conduct the present case study, it was decided to install the green coverage on the pre-existing roof of the building, made of hollow bricks, with a surface of approximately 80 m<sup>2</sup> .

As regards the weather conditions of the site, they were typical of the South of Italy, characterized by a temperate climate with warm summers and mild winters. Figures 2–4 show the trend of outdoor air temperature (T), relative humidity (RH) and solar radiation (IR), respectively, during the monitoring period of one year.

**Figure 2.** Trend of the outdoor air temperature during the monitoring period.

**Figure 3.** Trend of the outdoor air relative humidity during the monitoring period.

**Figure 4.** Trend of the solar radiation during the monitoring period.

#### *2.2. Description of the Analysed Green Roof Installed in the Experimental Site*

Going from the indoor to the outdoor sides of the building, as shown in Figure 5, the green roof compound is composed of the following layers: a root barrier (with a waterproofing membrane), a drainage layer (made of a polyethylene geo-net, hot-coupled with a non-woven geotextile with filtering functions), a water storage layer (constituted by cushions filled with expanded perlite), a filter fabric (composed by a geotextile felt, 100% polypropylene calendered), a growing medium (which is a mixture of peat, lapillus, pumice, zeolite and slow releasing fertilizers and is infesting weeds free) and the vegetation layer.

**Figure 5.** Sketch of the green-roof layers.

In order to also analyse the effects provided by different plant species and different thicknesses of the water storage layer on the green roof thermal behaviour, the roof was divided into six sectors, where three different Mediterranean autochthonous species (*Halimione Portulacoides*, *Rosmarinus O*ffi*cinalis Prostratus* and *Crithmum Maritimum*) and two different thicknesses of the water storage layer (10 cm for plots 1, 2, 3 – P1, P2, P3 – and 15 cm for plots 4, 5, 6 – P4, P5, P6 –) were used according to the scheme reported in Figure 6.

**Figure 6.** Scheme of the plant species planted in the different sectors.

A brief structural description of each layer is reported in Table 1.


**Table 1.** Description of the layers constituting the green roof plot.

Table 2 reports, instead, the main physical parameters characterizing each of the green roofs' six plots. The data listed in Table 2 are the same as those used in [69].


**Table 2.** Description of the layers constituting the green roof plot.

In addition to thermal conductivity, density and specific heat, which characterize the thermo-physical behaviour of the soil, some other properties typical of the specific plants, which have an important impact on the heat exchanges through the green-roof, were also considered.

Particularly, the "leaf reflectance" (dimensionless)**,** that is, the ratio of the incoming light which is reflected by a leaf, and the "leaf area index"—LAI (m<sup>2</sup> /m<sup>2</sup> )—defined as the one-sided green leaf area per unit of ground surface area. The latter, in particular, which has a great influence on the shading and transpiration effects, has a positive effect, especially during summer seasons; in fact, the higher the LAI, the higher the cooling reduction [50,69,70].

As for the vegetation characteristics, Figure 7 shows how the *Halimione Portulacoides* (both P1 and P6) and the *Rosmarinus O*ffi*cinalis Prostratus* (only P5) reached full coverage (100%) in less than 12 months while the *Rosmarinus O*ffi*cinalis Prostratus* in P2 achieved a maxim coverage of about 85% in the same period and the *Crithmum Maritimum* (both P3 and P4) did not accomplish more than 40%–60%, showing the difficulty of establishing it in the considered environment [52].

**Figure 7.** The six-plots green coverages system.

#### *2.3. Data Monitoring System Adopted*

The field experimental part of the proposed approach essentially consisted in a monitoring campaign of the ceiling temperature values of the building.

Since the monitoring of the temperatures profiles of the ceiling was aimed at checking the effects of the presence of the green roof on the indoor conditions, particularly in terms of thermal comfort levels, measures were performed in the center of the rooms' ceiling, far from thermal bridges and lights fixtures, as shown in Figure 8.

**Figure 8.** Layout of the building, green roof plots' arrangement and positions of probes for the temperature measurement.

Concerning the measuring method, temperatures were recorded with a sample rate of 10 minutes, during both the winter and the summer seasons, by means of insulated T type thermocouple probes.

The monitoring campaign lasted one year and started six months after the installation of the green roofs in order to have the green coverage well stabilized and to allow the testing of the acquisition system.

#### *2.4. Simulations Performed in the Study*

The modelling part of the proposed approach consisted in utilizing the very popular EnergyPlus simulation code to run the building's thermal calculations. For this purpose, different scenarios were implemented, specifically:


Regarding the HVAC schedule of Scenario #2 and Scenario #3, it was decided to use a simple on/off schedule, with the HVAC working between 7:00 and 17:00, considering the same power capacity as that used in Scenario #1.

The authors would like to underline here that Scenarios #1 and #2 are characterized by the presence of green roofs, while Scenario #3 refers to a standard roof. The difference consists in the fact that, while the simulation of Scenario #2 totally complies with the Energy Plus code (by utilizing its typical green roof simulation routine), Scenario #1 is modelled by forcing the Energy Plus code, that is, imposing the monitored indoor ceiling temperatures as boundary conditions. In this way, the model was driven with real data based on the presence of the experimental green roof, avoiding the actual simulation of the green roof element itself. Hence, this allowed us to compare the PMVs obtained from the Energy Plus green roof simulation tool (with its relative assumptions limits) with those calculated on the basis of the real monitored experimental data.

In order to assess the direct and indirect reduction of CO2 emissions, in this work, a value of 85 kgCO2/ha per year [71] was considered for the direct reduction of CO2 emission, based on the extension of green covering, while a value of 0.531 tCO2/MWh (the average emissions for the current electric Italian energy mix [34]) was used to estimate the indirect reduction of CO2 emission based on energy saving for climatization purposes and having set a COP value equal to 3 for the cooling season and 3.5 for the heating season.

#### **3. Results and Discussions**

In this section, the results relative to the monitoring campaign and to the energy performance simulations are reported.

#### *3.1. Monitored Data*

Table 3 shows the average temperatures measured on the ceiling of each room sited below the green roof's six plots, for summer (July) and winter (February) conditions, in periods during which the air-conditioning system was working.


**Table 3.** Monitoring results of the green roofs six plots.

The monitoring results point out a general tendency to attain lower temperatures when the green coverage is higher, i.e., P1 (*Halimione Portulacoides*). Indeed, this plot shows that ceiling temperatures were generally 1–3 ◦C lower with respect to the other plots in summer and 1–2 ◦C higher during winter, hence representing a benefit for both the summer and winter seasons.

Moreover, it must be noted that the LAI has a positive influence on the green-roof thermal behaviour; P1 and P6 have, in fact, higher LAI, unlike P3 and P4. In addition, another factor that could have influenced the obtained results is represented by the light colour of the plants' leaf surface, which enabless a higher amount of solar radiation to be reflected.

The results shown in Table 3 also highlight the influence of the different type of plants. In particular, *Halimione Portulacoides* (P1 and P6) reduces temperature peaks more consistently. Therefore, this type of plant seems to be more suitable for lowering the summer temperature values and increasing the winter temperature peaks.

Anyway, as reported in Table 3, it should be pointed out that during the summer season, a mean temperature of about 26.6 ◦C has been recorded by the thermocouples placed on the rooms' ceilings under P1 (*Halimione Portulacoides*), with maximum peaks of 27 ◦C, that lies within the suggested range for the indoor comfort in summer conditions [72]. The same cannot be stated for the other plots, where, even with a 100% green coverage, ceiling temperatures of about 27–28 ◦C were registered, with maximum peaks going beyond 30 ◦C.

The above-discussed outcomes demonstrate the importance of selecting a proper plant species during the green roof design phase.

Apart from these considerations, strictly related to the physical characteristics of the green roof, it is also necessary to take into consideration some aspects related to the building features that may have influenced the monitored ceiling temperatures, in particular:


#### *3.2. Outcomes of Energy Simulations*

Since the main aim of this work was to assess how the use of green roofs can affect the indoor comfort and the energy consumption of a building, it was decided to report, in the first part of this section, a comparison between the simulation results of Scenario #1 and Scenario #3. In particular, considering that such estimation is affected by the temperature changes during the actual HVAC system operating periods and in light of the detailed schedule utilized to run the simulations, it was chosen to divide the resulting data into two five-days periods representative of winter (Figure 9, on the left) and summer (Figure 9, on the right) conditions. Specifically, in Figure 9, the green lines represent an average of the values obtained for the six plots (Scenario #1), with its relative ranges of variation and the blue lines represent the standard case (Scenario #3); on the other hand, the red lines indicate the HVAC system start-up intervals.

**Figure 9.** Comparison between Scenario #1 (green lines) and Scenario #3 (blue lines) for winter (left) and summer (right) conditions.

Looking at Figure 9, it can be noted how after an initial start-up phase of the HVAC system, the presence of a green roof during the winter season does not seem to improve the indoor thermal conditions, while during the summer season, it brings a noticeable improvement of the indoor comfort levels. It must be underlined here that the amplitude of the variation range relative to the green roofs' temperature values is due to the fact that as reported in Section 2.2, the six plots are characterized by different features and therefore, describing parameters; in particular, the type of species and its relative coverage percentage have a strong influence on the monitored temperatures.

The temperature differences noticed also had an impact on the energy consumption. Over the representative five-day periods considered, in fact, a 18% increase for heating needs (220.05 kWh for the standard roof against 259.59 kWh for the green roof) and a 44% saving for cooling needs (189.38 kWh for the standard roof against 106.37 kWh for the green roof) were observed.

As mentioned earlier, after this first comparison, the authors wondered what results would have been obtained by simulating a green roof similar to the real one using the green roof configuration tool provided by EnergyPlus. The second part of this section shows, therefore, a comparison between the results obtained from the Scenario # 1 and Scenario # 2 simulations.

Similarly to the previously reported Figure 9, Figure 10 contains the obtained results for the two five-days periods representative of winter (on the left) and summer (on the right) conditions. In particular, the green lines (Scenario #1) and the red lines (HVAC system start-up intervals) are the same as shown in Figure 9, while the black lines represent an average of the values obtained for the six plots, with the relative ranges of variation, using the Scenario #2 EnergyPlus settings.

**Figure 10.** Comparison between Scenario #1 (green lines) and Scenario #2 (black lines) for winter (left) and summer (right) conditions.

By analysing Figure 10, it can be seen how, in the winter conditions, the green roof simulated according to Scenario #2 shows a very similar behaviour to that in Scenario #1, particularly during the air-conditioning working periods. In winter conditions, however, Scenario #2 allows to obtain higher temperatures than Scenario #1, corresponding to a further improvement in the indoor comfort levels. Furthermore, in summer conditions, contrarily to Scenario #1, Scenario #2 shows an evident very variable temperature trend between day and night, typical of a context highly influenced from solar radiation, which does not seem to reflect reality.

As for the fact that the changes of the Scenario #2 temperatures range are much narrower than those in Scenario #1, it must be observed that this is most likely due to the fact that, as previously highlighted, the model used by EnergyPlus does not allow to set all the parameters of the green roof freely but imposes some constraints to their numerical values. Due to this reason, in fact, the Scenario #2 results show no differences relating to the two different thicknesses of water storage used for each species, but only some small differences between the different species.

As for the energy consumption of Scenario #2, over the representative five-day periods considered, a 4% saving for heating needs (220.05 kWh for the standard roof against 212.34 kWh for the green roof) and a 41% saving for cooling needs (189.38 kWh for the standard roof against 112.54 kWh for the green roof) were observed.

Finally, in the last part of this section, it was decided to report a rough estimate, on a monthly basis, relative to both aspects of indoor comfort improvement and energy consumption savings. In this regard, it was chosen to compare the results deriving from the simulations of Scenarios #2 and #3. This choice was suggested due to the fact that, given the intended use of the building (i.e., research laboratory) and its real use (i.e., infrequent), especially in terms of the HVAC system, it was assumed that the results of Scenario #1 were not considered actually representative for a long-term estimate.

Once again, in order to make the results visually more easily readable, an average of the results obtained for the green roof six plots was used to display the results of Scenario #2.

To compare Scenario #2 with Scenario #3 in relation to indoor comfort levels, it was decided to report, in Figure 11, a graph of the monthly temperatures. The graph indeed allows to show the deviations of the average values of the green roof ceiling temperatures (∆tav,ceiling) compared to those of the standard roof, where the black bars represent the range of deviation from the average values.

**Figure 11.** Deviations of the average values of the green roof ceiling temperatures compared to those of the standard roof.

Figure 11 highlights the positive effects due to the presence of the green roof, which, with respect to the standard roof, allows maintaining higher ceiling temperatures in winter and lower ceiling temperatures during summer.

Other than the temperature, another important indicator when assessing the indoor comfort levels is represented by the PMV (Predicted Mean Vote). For this reason, it was also decided to report, in Figure 12, a comparison between the monthly PMV average and peak values of Scenario #2 (GR) and Scenario #3 (ST).

**Figure 12.** Monthly PMV average and peak values for Scenario #2 (GR) and Scenario #3 (ST).

By looking at the differences in the obtained PMV average values (Figure 12) with and without the presence of the green roof, especially those relative to July and August, it arises that the presence of the green roof reduces PMV average values from more than 0.7 to approximately 0.5. Hence, accordingly to the standard currently in force for the design of the indoor environment, i.e., EN 16798-1:2019 [72], the presence of the green roof contributes to shift the indoor thermal environmental conditions from Category III (acceptable, moderate level of expectations) to Category II (normal level of expectation). In other words, the presence of the green roof contributes to bring the building within comfort conditions (PMV = 0.5), starting from a slight warm condition (PMV = 0.7).

Moreover, by analysing Figure 12, it can also be seen how, although a general positive effect due to the presence of the green roof is evident, some critical issues emerged in the months of April and October (transition months), for which the standard roof seems to perform better than the green roof. This condition, which needs to be better investigated, is probably due to the additional thermal inertia that the presence of the green roof brings to the structure: this slows down the response of the green roof compound to the changes of climatic conditions occurring in the transition periods of the end of spring (April) and the beginning of winter (October).

For the sake of completeness, it was decided to report, in Figure 13, an annual plot where the average daily external temperatures (Outdoor) are compared with those of the ceiling in the presence of the green roof (GR\_mean) and with those relative to the standard roof (ST).

**Figure 13.** Annual average daily temperatures (i.e. Toutdoor air, Tceiling with GR and Tceiling without GR) trends.

Regarding the energy consumptions for heating and cooling needs, these are summarized in Table 4 by reporting the absolute values (kWh) obtained for the standard roof scenario and the correlated average percentage deviations (including the respective variation ranges) relative to the achievable savings due to the green roof presence.


**Table 4.** Monitoring results of the green roofs six plots.

Specifically, for each month, the amount energy consumed for both heating and cooling with and without the presence of the green roof is listed. The use of the bold character is intended to show more easily the actual HVAC working periods, that is December–March for winter conditions and June–September for summer conditions. Therefore, in Table 4, data related to months when the HVAC is working have been highlighted using the colour black. The results confirm the advantage of using green roofs as a solution capable of achieving valuable energy savings.

Moreover, the mean indirect reduction of CO<sup>2</sup> emissions due to the green roof installation was 145.6 ± 13.8 tCO2/year, while for the direct reduction, a value of only 56 gCO2/year was observed.

Finally, some further considerations need reporting. When an energy restoration of a roof is in context, such as the one considered in this work, it is easier and safer to add a vegetated coverage to an existing roof than making it larger. That is why here, it was chosen to exclude a theoretical comparison between a green roof and a hypothetical alternative high massive one and to limit the analysis to a specific comparison between the behaviour of the pre-existent standard roof equipping the building and the improvements brought by the installation of the green coverage. For this purpose, the thermo-physical characteristics of the roof are those typical of the building habit of the considered geographical area.

In addition, the benefits of the presence of a green roof cannot be simply evaluated in terms of thermal insulation, since it generates other positive effects regarding the evaporative phenomena and the change of the albedo of the roof. Therefore, it was decided to exclude a comparison with insulated roofs too, in accordance with existing literature studies—such as that by Niachou et al. [73]—demonstrating that the presence of green roof on an insulated roof is practically irrelevant.

#### **4. Conclusions**

The capability of green roofs in reducing the electric energy needs for the climatization of buildings and their environmental benefits has been extensively demonstrated in the literature.

The idea behind the presented work derived from considerations regarding the possible influence of green roofs on the indoor thermal comfort levels; that is, instead, an aspect often overlooked when estimating the performance of such building components. Therefore, the analysis methodology approach utilized (partly modelling and partly experimental) was implemented in light of such considerations.

Accurate knowledge of the internal temperatures of the ceiling is indeed an important prerequisite for establishing both achievable indoor comfort conditions and the energy demand for the air conditioning of the building itself. A lowering of the indoor temperatures in summer, and a rise during winter, lead, in fact to an improvement in the comfort levels for the occupants, and consequently, a saving on the use of the HVAC system, which, in turn, translates into a reduction of polluting emissions.

For this purpose, the comparison between Scenario #1 (where the monitored ceiling temperatures were used as boundary conditions) and Scenario #3 (standard roof case), and that between Scenario #2 (where the green roof configuration provided by EnergyPlus was utilized) and Scenario #3, allowed to prove how the presence of the experimental green roof on the monitored building improved the indoor comfort levels during summer by moderating the ceiling temperatures (Figures 9–11 and 13), despite some worsening during winter periods seeming to occur. Moreover, the temperature differences noticed had also a positive impact on the building energy consumption and the CO<sup>2</sup> emissions.

On the other hand, by comparing Scenario #1 with Scenario #2, it was possible to highlight the possible limits of the code's ability in adequately simulating the green roof behaviour. These limits are mainly represented by the lack of flexibility that EnergyPlus allows in the setting-process of some of the green roof physical parameters and in the way in which the code takes into account the solar radiation components.

Another aspect which should be better investigated, regards the PMV results (Figure 12). In fact, even though, also, in this case, a general positive effect due to the presence of the green roof can be seen, some criticalities emerged during some transition periods, for which the standard roof seems to perform slightly better than the green roof.

In conclusion, the proposed analysis made it possible to highlight how it is possible to assess the impact that green roofs have specifically on the indoor comfort levels, other than on the energy consumption.

Furthermore, the availability of field data put into evidence the importance of adequate simulation tools to facilitate the green roof design and assessment processes.

**Author Contributions:** All authors contributed equally to the design, experimental analyses and editing of this research. Conceptualization, L.C., M.L.G., G.P., G.R., G.S. (Gianluca Scaccianoce), G.S. (Giancarlo Sorrentino) and S.A.; Data curation, L.C., M.L.G., G.P., G.R., G.S. (Gianluca Scaccianoce), G.S. (Giancarlo Sorrentino) and S.A.; Methodology, L.C., M.L.G., G.P., G.R., G.S. (Gianluca Scaccianoce), G.S. (Giancarlo Sorrentino) and S.A.; Writing—original draft, L.C., M.L.G., G.P., G.R., G.S. (Gianluca Scaccianoce), G.S. (Giancarlo Sorrentino) and S.A.; Writing—review and editing, L.C., M.L.G., G.P., G.R., G.S. (Gianluca Scaccianoce), G.S. (Giancarlo Sorrentino) and S.A. 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.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Towards Nearly Zero Energy and Environmentally Sustainable Agritourisms: The E**ff**ectiveness of the Application of the European Ecolabel Brand**

**Laura Cirrincione 1,\* , Maria La Gennusa <sup>1</sup> , Giorgia Peri <sup>1</sup> , Gianfranco Rizzo <sup>1</sup> and Gianluca Scaccianoce 1,2**


Received: 28 July 2020; Accepted: 17 August 2020; Published: 19 August 2020

**Abstract:** Tourism represents an important economic driver in Italy, being responsible for approximately 13.2% of the total GDP (a value higher than the reference European average) and for nearly 10% of the regional GDP. Among the touristic sectors, the agritourist ones show a persistent growth, experiencing in 2019 a 6.7 point percentage improvement compared to the 2017 figures. Given this situation, the transition towards a low-carbon path, affecting the building sector for some time, should also involve agritourist buildings, through the release of EU directives, member state laws, and technical rules. On the other hand, agritourism sites could be awarded the Community EU Ecolabel. Unfortunately, awarding the EU environmental excellence brand implies the availability of several data on building energy behavior that should then be managed by complex evaluation tools. To overcome this issue, the use of the simplified ARERA (Italian Regulatory Authority for Energy Networks and Environment) technical datasheets, issued to assess environmental improvements consequent to energy efficiency interventions in the urban residential building stock, is proposed. The application of this tool totally avoids using building computer-based simulation models, thus facilitating the preparation of the EU Ecolabel request documentation by agritourism owners. Being awarded the Community EU Ecolabel also implies approaching a net zero energy condition because of a lower energy consumption and a minor recourse to fossil fuels. For this purpose, an application of an easy graphical method, previously developed for residential and commercial buildings, which visually represents improvements achievable by a given agritourism when implementing energy efficiency measures, is presented.

**Keywords:** building energy efficiency; European environmental brands; tourism sector; agritourism; nearly zero energy buildings (nZEB)

#### **1. Introduction**

Tourism, and the activities connected to it, represent an important sector of the economic system. According to recent statistics the tourism industry represents about 10% of total global gross domestic product (GDP) and 7% of global trade [1,2], accounting for approximately 11% of the world's employment, with an expected positive economic growth trend [3,4]. Tourism constitutes a significant contributor to energy consumption, both at a global and European scale [5–7], which translates to a significant impact on the environment and ecosystem; it is in fact responsible for about 5% of the global CO<sup>2</sup> emitted by human activities [1,8].

Accommodation (thus the building), in particular, is the third energy consuming item (after travel and transport), much of which is consumed in space heating or air conditioning (up to 50% in some cases), followed by hot water, and cooking [9,10]. Moreover, a study conducted by the World Tourism Organization and the United Nations Environment Programme [11] estimates that accommodation generates 21% of tourism's total greenhouse gas (GHG) emissions. Accordingly, the number of papers analyzing tourism significance, in terms of energy consumption [12,13] and impacts on emissions [14,15], has been increasing lately.

Consequently, in recent years much attention has been paid to the concept of sustainable tourism, which in accordance with the United Nations Environment Programme (UNEP) and the United Nations World Tourism Organization (UNTWO) is defined as "development of tourism activities with a suitable balance between the dimensions of environmental, economic, and sociocultural aspects to guarantee its long-term sustainability". Hence, the challenge of sustainable tourism is to mitigate its negative impacts, consisting mainly in: (i) high energy consumption, (ii) increasing GHG emissions in the atmosphere, and (iii) the contribution to climate change [16].

Therefore, taking into consideration global [17,18], European [19–22], and national [23] policies, the UNWTO recommended three central actions on which the tourism sector should concentrate in order to contribute in achieving a more sustainable development [1,24], which are resource efficiency, environmental protection, and climate change (linked to sustainable development goals (SDGs) 6, 7, 8, 11, 12, 13, 14, and 15) [25].

At the European scale, the European Commission set the basis for the best environmental management practice in the tourism sector in accordance with Article 46 of the Eco-Management and Audit Scheme (EMAS) regulation [26,27]. Furthermore, by means of the "Guide on EU funding for the tourism sector 2014–2020" [28] the EU states that effective governance, policies, frameworks, and tools need to be implemented in order to properly guide and support (also from an economic point of view) the development and promotion of sustainable tourism practices.

Tools like these are indeed important because they encourage the owners and/or managers of the accommodation facilities to use practices and systems that allow both energy savings and pollutant emissions, by favoring the visibility of these structures in terms of environmental sustainability, which represent an added value, given that tourists are becoming increasingly more attentive to this issue.

In this regard one of the first initiatives undertaken by the European Community has been the releasing of the EU Ecolabel for tourist accommodation services [29], created to improve the environmental performance of hotels, campsites, hostels, agritourisms, holiday homes, and bed & breakfasts, by providing efficient guidelines on the action to be implemented in order to lower their environmental impact; and which still remains one of the most implemented initiatives.

The promotion of sustainable tourism is also the basis of the nearly zero-energy hotels (neZEH) project, launched by the Intelligent Energy Europe Programme of the European Commission, with the intent of supporting European hotels in complying with the nZEB (nearly-Zero Energy Buildings) regulations [1]. On this subject, various studies have been conducted aimed at analyzing the achievable energy saving measures [30–32] and proposing suitable strategies and policies to be adopted [33,34].

Looking at the national scenario, the tourism issue is particularly relevant, considering that 16.5% of EU accommodation facilities are located in Italy [35], and since in the last two years Italy resulted to be amongst the top five most visited European tourist destinations (for accommodation in hospitality facilities), with a 13.4% share of the total of the EU-28 [36,37].

According to some recent statistics, the Italian tourism sector represents 13.2% of the national GDP (for a total contribution of around 230 billion euros), higher than both the world and European figures (which stand at around 10%). The economic impact of tourism is significantly reflected in the job market, accounting for 14.9% of the country's total employment [38]. Tourism is in fact one of the fastest growing industrys in Italy, and both public and private business organizations are strongly interested in its economic and environmental impact, both at national and regional level [39,40].

Thus, from the collaboration between such organizations and the national and regional governments, different initiatives have been undertaken from an environmental sustainability point of view in recent years. These include the creation of a set of national environmental quality certifications (besides the previously cited EU Ecolabel), including the "Green Key" [41], "Bandiera Blu" [42], and "Spighe Verdi" [43], born from the collaboration between the Italian Foundation for Environmental Education—FEE Italia (whose actions are supported by ONU, UNEP, UNWTO, and UNESCO) and national authorities dealing with environmental policies [44].

Furthermore, other economic initiatives have been implemented to encourage the use of sustainable energy solutions through financial incentives. The "Tax Credit Alberghi—Bonus alberghi e agriturismi" (bonus for hotels and agritourism), a tax facility that encourages various upgrading activities, including those aimed at improving energy efficiency, has recently been introduced, specifically for accommodation facilities [45].

Agritourism, or rural tourism, has been promoted as a practice able to encourage the use of green practices, making farms sustainable and also maintaining the local historical and natural settings [8,46].

Thanks to this, according to recent statistics in Italy, the agritourism sector continues to record a growing trend, both in the number of structures, and in the presence of customers and its economic value. Agriculture economic reports make it possible to measure the economic dimension of the agritourism sector, which is equal to 1.36 billion euros, up 6.7% compared to the previous year. In particular, 60% of agritourisms are located in the regions of central and southern Italy, where Sicily prevails with more than 600 farms [47].

The growing interest in the agritourism sector is also reflected in the academic world, where studies concerning both the economic and social benefits of various tourist activities in the rural area, including agritourism [48], and the environmental performance of agritourism companies in terms of energy performance [49,50], can be found.

In the present work we verified whether the simplified ARERA (Italian Regulatory Authority for Energy Networks and Environment) technical datasheets [51], issued for the urban residential building stock, can be easily applied to estimate the increase in energy efficiency (or the corresponding decrease in the release of polluting substances) consequent to the adoption of some improvements to a building or plant, planned for the issuance of the EU Ecolabel brand for accommodation facilities [29]. The convenience in the use of these technical datasheets lies in the fact that they allow the estimation of the energy demand reductions without necessarily going through the building simulation. For this purpose, a case study has been conducted to estimate what advantages agritourism owners could gain in adopting a well-known brand such as the EU Ecolabel [29], with particular reference to the actions aimed at saving energy and reducing emissions of pollutants, from the perspective of a possible "nearly Zero Energy Agritourism (nZEA)", in parallel with the previously cited nZEB and neZEH projects.

The idea at the base of this work stems from the numerical consistency of agritourisms in Sicily [47] and their conceivable growth trend, which is a consequence of the increased interest in the rural landscape of the territory and in the products of the land that are strongly orienting tourism, directing it not only towards the urban context. The adoption of an environmental certificate like the EU Ecolabel [29] can therefore represent an advantage both for agritourism owners and for the entire territory.

Furthermore, the owners of agritourism in Sicily can apply for subsidized loans and financial funding [28,52] in the regional area and beyond. However, such requests must be supported by information concerning the consumption and energy efficiency of the agritourism and, in line with the new European directives on sustainability [17–19], by information on the environmental performance of the buildings themselves (premises).

Normally this information is of a complex nature and tends to imply the use of sophisticated simulation models, the use of which is not always the prerogative of (or available to) the managers of the holiday farms. The same problem can be found by analyzing the work of the decision makers who have to assess the adequacy of the requests for funding.

Essentially, the availability of simple but reliable tools for evaluating these premises is of paramount importance for the orientation of this important tourism sector towards a sustainable path.

Hence, as previously mentioned, in order to provide a contribution to this important issue we assessed the reliability of a scheme of simple computational methods provided by the Italian Regulatory Authority for Energy Networks and Environment—ARERA [51], specifically for the residential and tertiary building stock. The advantage in the use of this computational scheme lies indeed in the fact that it is based on excel spreadsheets (technical datasheets) which, as already mentioned, allow the estimation of the energy demand reductions without the need of simulating the building behavior.

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

The proposed methodology aims at considering together in an easy and accessible way two aspects of the sustainability, which are energy efficiency and environmental safety, in order to help agritourism owners, and/or managers, to make decisions that are more favorable to them and consistent with the European policies in force. Specifically, according to the presented approach, the selection of energy efficiency interventions is based on a combination of the ARERA technical datasheets and the EU Ecolabel criteria, hence taking into account the environmental sustainability aspects, and also in view of achieving a possible nearly zero energy condition (nZEA). Therefore, two Sicilian agritourisms have been selected to show how the application of the proposed methodology actually works.

The considered approach can also be seen as a simple diagnosis method aimed at facilitating the social appropriation of knowledge and technology, so that the owners of agritourism facilities can confidently check their level of eco-efficiency. Moreover, the method can be utilized in order to choose between addressing actions concerning the energy performance of the structure or interventions regarding the installation of new (renewable) energy plants.

#### *2.1. Agritourism Definition*

The Italian national legislation [53], and the regional Sicilian one [54,55], define as 'agritourism' activities, those reception and hospitality activities exercised by agricultural entrepreneurs, through the use of their own company connected with the activities of cultivation of the land, forestry, and animal breeding. Thus, agritourism activities include:


#### *2.2. The ARERA Technical Datasheets*

In the present work, the use of ARERA technical datasheets [51] was not an arbitrary (random) choice, but it was decided to turn to these methods since, although simplified, they constitute an official reference at the Italian national level.

The Italian Regulatory Authority for Energy Networks and Environment—ARERA is indeed an independent body, established with the task of protecting consumers' interests and promoting competition, efficiency, and the spread of services, and having adequate quality levels, through regulation and control activities. The action of ARERA concerns the sectors of electricity and natural gas [56], water services [57], district heating and district cooling [58], and the waste cycle [59].

One of the main tasks of ARERA is to promote the rational use of energy, with particular reference to the promotion and diffusion of end use energy efficiency and/or energy saving actions, and the adoption of measures for sustainable development. Among the feasible actions, there are both active measures, which involve the installation of high efficiency equipment, or the insertion of regulation devices for a more efficient use of energy, and passive interventions such as the modification of buildings' envelope in order to reduce losses. modification of buildings' envelope in order to reduce losses.

In this regard, the technical datasheets proposed by ARERA establish the guidelines for the preparation, execution, and final evaluation of specific actions, aimed at increasing energy efficiency (or promoting energy saving), providing reduced rates of primary energy consumption actually achieved (expressed in toe—Tons of oil equivalent), and also for the purpose of issuing energy efficiency certificates. Table A1 in Appendix A reports a comprehensive list of the current standardized and analytical ARERA technical datasheets. —

#### *2.3. The EU Ecolabel Brand*

Established in 1992 (by Regulation n. 880/92 [60], now disciplined by Regulation (EC) n. 66/2010 [61] in force in the EU-28) and recognized across Europe and worldwide (Figure 1), the EU Ecolabel is a voluntary environmental performance certificate that is awarded to products and services meeting high environmental standards. The EU Ecolabel encourages companies to develop products and provide services that consume less energy, and generate less waste and CO<sup>2</sup> emissions. As of March 2019, an increase by 88%, with respect to 2016, of the number of EU Ecolabelled products/services has been registered. Leading countries for number of products/services are: Spain, Italy, Germany, Belgium, and France [62].

**Figure 1.** Official EU Ecolabel logo [62].

The EU Ecolabel provides exigent criteria, and relative guidelines, depending on the type of product and/or service, in order to reduce their overall environmental impact. Such criteria are established at a European scale with a wide participation of interested parties, including both public authorities, and consumer and environmental associations [63].

In particular, the EU Ecolabel for tourist accommodation services [29] was created specifically for hotels, campsites, hostels, agritourisms, holiday homes, and bed & breakfasts, in order to improve their environmental performance, by providing a set of criteria on the action to be implemented in order to lower their impact. Such criteria are divided into mandatory and optional, and focus on the five categories; general management, energy, water, waste and wastewater, and other, as shown in Table A2 in Appendix A.

within the product group "tourist accommodation" according to the legal obligations of the country In order to be awarded the EU Ecolabel a tourist accommodation service, other than falling within the product group "tourist accommodation" according to the legal obligations of the country in which the accommodation is located, must comply with all the mandatory criteria (if applicable), and receive at least twenty points under the optional criteria [29].

as a way of legitimizing the accommodation's An added value, in terms of visibility, for tourist accommodation owners lies in the fact that the EU Ecolabel is recognized by the majority of travelers as a way of legitimizing the accommodation's claims that it is making real efforts to reduce its impact on the environment in its operational activities.

#### *2.4. Analysis Methodology*

#### 2.4.1. Merging the ARERA Technical Datasheets and the EU Ecolabel Criteria

As mentioned in the introduction section, the aim of the present work is to verify whether the simplified ARERA technical datasheets can be applied to estimate the increase in energy efficiency consequent to the adoption of some actions planned for the issuance of a EU Ecolabel for Tourist Accommodation Services, without necessarily going through the building simulation. Obviously, an increase in the energy efficiency implies a corresponding decrease in the release of polluting substances.

Therefore, since this work is mainly focused on the energy criteria, starting from the assumption that the considered agritourism meets all the mandatory criteria, it was decided to analyze possible "environmental action packages", consisting of different combinations of the actions established by the optional energy criteria, which are better suited to a scenario such as agritourism, and which allow the obtaining of the weight of the energy category on the twenty points minimum limit set by the regulation, which corresponds to 7.34 points.

To this purpose, only the ARERA datasheets regarding the actions related to the improvement of the structure energy efficiency that could be transferred and applied to agritourism structures, according to the EU Ecolabel for Tourist Accommodation energy criteria, have been considered, as reported in Table 1, where the correspondent energy consumption categories have also been reported.


**Table 1.** Correspondence between the Italian Regulatory Authority for Energy Networks and Environment (ARERA) datasheets and the EU Ecolabel energy criteria.

<sup>1</sup> HVAC—Heating, Ventilation, and Air Conditioning; DHW—Domestic Hot Water; RES—Renewable Energy Source.

As can be seen in Table 1, under the dotted line, the two ARERA datasheets 6T and 20T, additional to those that can be associated to the EU Ecolabel, have also been taken into account. In fact, even though these two intervention typologies are not foreseen by the current Ecolabel scheme, they represent actions that can actually be applied to an agritourism structure in the perspective of a possible "nearly Zero Energy Agritourism (nZEA)" as a parallel with the well-known nZEB concept; and also, in view of a possible future improvement of the Ecolabel scheme.

The selected datasheets (Table 1), have been then put into the form of appropriate excel spreadsheets.

The equations relative to the ARERA calculation procedures, for each considered technical datasheet, are given in Appendix B.

One aspect that must be highlighted here regards the fact that while datasheets 5, 6, 15T, 19T, 20T, and 27T enable obtaining savings of consumed energy (energy saving measures, ESM), datasheets 7 and 8T allow, instead, the production of energy from renewable sources (renewable energy sources, RES).

#### 2.4.2. Methodology Application Feasibility

With the aim of assessing the potential energy savings, with reference to a real context, it was decided to select two agritourisms situated in the Sicilian province of Palermo, considered as representative of the whole regional agritourism context regarding the size, the provided services, and more importantly for the purpose of the proposed methodology object of the present work, in terms of energy consumption. Apart from these physical and energy features, both agritourisms were selected thanks to their wide offer of services, which are representative of these kind of farms, and due to the fact that they operate in the two climatic zones where agritourisms are mainly sited in

Sicily. Specifically, the two agritourism are *Villa Dafne*, sited in Alia and belonging to climatic zone D, and *Bergi*, located in Castelbuono and classified as climatic zone C. Both agritourisms fall into solar belt 3. Table 2 describes the main general characteristics of the two structures.


**Table 2.** General characteristics of the two selected agritourisms.

In Table 3 is reported the information relative to the energy characteristics of interest for the conducted study, which were obtained by on field surveys and interviews with the owners of the two businesses, thanks to which it was possible to reconstruct the energy consumption relative to an entire year of operation of the structures. In particular, the data regarding the energy consumption were distributed between four main categories and accordingly broken down into percentages, and corresponding toe/year, also with reference to the corresponding energy sources. As for the energy sources' average costs, the following values were used:



**Table 3.** Energy sources and energy consumption breakdown for the two selected agritourisms.

Subsequently, it was hence possible to obtain the achievable energy savings (AES), in terms of percentage of electricity consumption covered by the datasheets proposed interventions on an annual basis, by comparing the values obtained from Equations (A1) to (A8), and the total energy consumption (Table 3), by means of the following equation:

$$AES\_i = \frac{R\_i}{Tot.\\_cons\_j} \text{ [\%]} \tag{1}$$

where:


Regarding the pollutant emissions, an assessment of the CO<sup>2</sup> emissions' reduction was conducted assuming for the considered climatic context an emission factor equal to 2.30 tCO2eq/toe for the electrical supply [64,65], while for natural gas and diesel oil an emission factor of 3.08 tCO2eq/toe and 2.34 tCO2eq/toe, respectively [66].

In order to single out the most convenient aforementioned "environmental actions packages", an economic estimation relative to the interventions suggested by the ARERA technical datasheets was also performed. To this purpose, the information relative to the costs of supply and installation for the materials, used to calculate the proposed interventions costs, were obtained from the current regional price list [67] and from local market surveys, as reported in Table 4.


**Table 4.** ARERA technical datasheets proposed interventions costs.

Successively, the economic savings, in terms of saved €/year, were obtained by multiplying the energy savings with the energy sources' average costs, according to the considered categories breakdown (Table 1). Furthermore, in order to select the optimal "environmental actions packages", for these the pay-back periods (not discounted) were also calculated and expressed in years.

#### **3. Results**

In this section the outcomes of the application of analysis methodology are reported.

Regarding the input parameters used in the equations relative to the ARERA calculation procedures, for each considered technical datasheets, these are given in Table A3 in Appendix B.

The following Figures 2 and 3 show the achievable energy savings (AES) on the total annual consumption, relative to the application of the intervention proposed by each considered ARERA datasheet, to the two agritourisms. On the right side of the graphs, the EU Ecolabel points corresponding to each datasheet are also reported.

*Villa Dafne* agritourism' achievable energy savings for each datasheet

**Figure 2.** Achievable energy savings (AES), on an annual basis, and EU Ecolabel points relative to each considered ARERA datasheet for *Villa Dafne* agritourism.

*Bergi* agritourism' achievable energy savings for each datasheet

**Figure 3.** Achievable energy savings (AES), on an annual basis, and EU Ecolabel points relative to each considered ARERA datasheet for *Bergi* agritourism.

Figures 2 and 3 show how, amongst the ARERA proposed interventions enabling the obtaining of an EU Ecolabel score, the implementation of a photovoltaic (PV) system (datasheet 7) would be the one allowing the gain of a greater advantage in terms of energy consumption. Concerning the savings related to the domestic hot water (DHW) category, by comparing datasheets 8T and 27T, which are alternatives to each other, it can be observed how 8T would be the more convenient choice. Regarding the heating, ventilation, and air conditioning (HVAC) category, the savings achievable through the application of datasheets 15T and 19T should, instead, be considered jointly (15T + 19T) as they can be attributed to the same improvement intervention. As for datasheet 5, although it represents one of the easiest measures to implement, it does not seem to bring the benefits that would have been expected.

With respect to datasheets 6 and 20T, also in this case the consideration that the attainable benefits must be considered together (6 + 20T) is valid. As already explained these two datasheets fall out of the Ecolabel scoring scheme, nevertheless they represent the second-best intervention that allows the highest energy savings, after datasheet 7.

The overall obtained results for the two considered agritourisms are reported in Tables 5 and 6.


**Table 5.** Proposed interventions costs and environmental benefits for *Villa Dafne* agritourism.

**Table 6.** Proposed interventions costs and environmental benefits for *Bergi* agritourism.


As can be observed, according to what has been previously pointed out, a single intervention cost was given to datasheets 15T and 19T as the proposed intervention corresponds to the same type of system, i.e., the same system allows operation for both heating and cooling. The same consideration can be made for datasheets 6 and 20T in relation to the insulation of the building.

Looking at the economic savings column it can be noticed how, from this point of view greater advantages can be associated with datasheets 7 and 8T, followed by 15T + 19T, 27T and lastly 5. Considering the whole set of interventions, instead, the (6 + 20T) option would also result second in this case.

Referring to CO<sup>2</sup> emissions reduction, the obtained results are obviously in line with what was seen beforehand (Figures 2 and 3) and commented on with the energy savings.

#### **4. Discussion**

The application of the ARERA data sheets to agritourism raises a question concerning the suitability of these simplified forms to the energy performances of agritourisms sites, being originally developed for residential and commercial buildings.

On the other hand, a possible improvement of the energy features of an agritourism, due to the actions referred to in the ARERA datasheets, should be evaluated on the base of its effectiveness in addressing a given site, towards a nearly-zero energy path, as required by the current international standards [20,21,28].

Both issues are briefly discussed in the following.

#### *4.1. E*ff*ectiveness of the Proposed Actions*

The obtained results could seem not too encouraging in terms of energy savings. In fact, the reduction of the energy demand following the proposed actions accounts for about one third of the annual energy consumption for both the considered agritourisms. However, this is not surprising; the fact that the ARERA technical datasheets proposed interventions have been designed for the residential sector, in fact, place some limits on their application in a wider context, such as the agritourism one. Specifically, the limitations set on the reference physical units (UFRs) sizes might have made the outcomes much lower than the actually achievable results.

For instance, concerning datasheet 7 a maximum *kW<sup>p</sup>* of 20 kW is reductive for an agritourism, which could employ PV better having wide areas available to install such systems. Supporting this observation, during the survey of the agritourisms, it arose that both currently have a 100-kW PV undergoing design phase. In this context it would be more sensible to impose a limit on the maximum percentage of yearly energy consumption to be covered with the proposed intervention.

The latter consideration also applies to datasheet 8T.

Regarding, instead, datasheet 15T the application problem is mainly related to the residential standard apartment size (80–90 m<sup>2</sup> ), which is difficult to translate into an agritourism setting. In the conducted analysis, for instance, in order to comply with such a parameter, three to four rooms were grouped and assumed equal to 1.5 standard apartments, but it could be a questionable criterion.

As for datasheets 19T, it would be more reasonable to install a centralized system rather than considering the replacement of the single air conditioning units (the same goes for the heat pumps proposed by datasheets 15T).

Nevertheless, since one of the aims of this work was that of singling the most convenient EU Ecolabel "environmental actions packages", based on the comparison of the results reported in Tables 5 and 6 it was decided to tentatively choose three alternative options, both for *Villa Dafne* (*VD-n*) and *Bergi* (*B-n*), as follows:


Table 7 summarizes the obtained results relative to the selected "environmental actions packages". By analyzing the data reported in Table 7 it was, therefore decided to consider as optimal options *VD-2* for *Villa Dafne* and *B-2* for *Bergi*. These two options allow, in fact, the obtaining of greater economic and energy savings and, correspondingly, higher CO<sup>2</sup> emissions reductions. Moreover, they are characterized by the lower pay back periods.

It must be observed that the availability of effective and reliable methods for evaluating the energy actions involving agritourism is of paramount importance for suitable planning of this important sector. Therefore, the ARERA technical data sheets should be properly reconsidered in order to render them more complicit with the energy features of agritourism buildings and dwellings.


**Table 7.** Summarized results for the two agritourisms.

#### *4.2. Towards a Nearly Zero Energy Agritourism*

As already mentioned, another intention of this work concerned the possibility of applying some actions to the agritourism structures, additional to those envisioned by the EU Ecolabel, in order to move towards a potential "nearly Zero Energy Agritourism (nZEA)". For this purpose, the results relative to ARERA datasheets 6 and 20T were added to the selected optimal options *VD-2* and *B-2* for *Villa Dafne* and *Bergi*, respectively; the outcomes of such combinations are reported in Figures 4 and 5. move towards a potential "nearly Zero Energy Agritourism (nZEA)".

**Figure 4.** Path towards a nearly zero energy condition (nZEA) for *Villa Dafne*.

**Figure 5.** Path towards a nearly zero energy condition (nZEA) for *Bergi*.

– According to such approaches and visual representations, already used in the literature [68–70], the nearly zero energy condition (nZEA) is reached when the energy demand (reported on the x axis) is completely covered by the self-energy supply from renewable sources (reported on the y axis).

— — Therefore, the effectiveness of the ARERA proposed interventions in moving agritourism towards a sustainable path, nZEA, is given as a simple summation of the effects provided by the energy saving measures—EESM (datasheets 6, 15T, 19T, and 20T) and those attributable to the renewable energy sources—ERES (datasheets 7 and 8T). E<sup>A</sup> represents, instead, the current energy consumption and, D<sup>a</sup> and D<sup>p</sup> the current (ante operam) and achievable (post operam) minimum distances (hence the perpendicularity) from the nZEA condition, respectively.

For the sake of simplicity, this assumption does not take into account the synergetic effects that are likely induced by the contemporary adoption of different energy actions on a given agritourism site.

Consequently, the results reported in Figures 4 and 5 show that the selected combinations of interventions allow an improvement of 59% for *Villa Dafne* and a 62% for *Bergi*, in terms of approaching the nZEA condition with respect to the current conditions.

and ranking the "environmental actions packages" applicable to agritourisms Regardless of the obtained results, the proposed methodology can be seen as a simplified scheme for analyzing and ranking the "environmental actions packages" applicable to agritourisms, and could be usefully adopted by local administrations to define the impact of different scenarios in order to better define environmental policies concerning the agritourism sector.

The proposed assessment/estimation methodology could, therefore, also represent important information for the design of the rural tourism sector, and of the/a regional energy plan by stakeholders and decision makers [71,72].

#### *4.3. On the Correspondence between the EU Ecolabel Criteria and the ARERA Technical Datasheets*

The application of the ARERA technical datasheets and the EU Ecolabel criteria to two different agritourisms in Sicily (here considered representative of the whole regional agritourism context) enabled us to better understand the level of compliance between two such schemes. The level of correspondence cannot totally match, since the ARERA methodology has been designed specifically for residential buildings and, on the other hand, the EU Ecolabel for Tourist Accommodation Services applies expressly to tourism facilities, with features that could not be perfectly applicable to agritourisms. These latter, in fact, are generally characterized by the presence of cultivated soils and the production of agrifarm foods and products.

Nevertheless, the comparison exerted on the two sites has shown that some useful correspondences can be assessed. In fact, by means of the combination of the ARERA technical datasheets and the EU Ecolabel energy optional criteria, it is likely possible to identify some "environmental actions packages", suitable to the agritourism context. Such packages are allowed to obtain a 7.34 points minimum limit for the energy category, set by European regulation. In particular, it emerged that the combination of datasheets 5 and 7 (options *VD-1* and *B-1*) allowed obtaining 7.5 EU Ecolabel points, while by adding datasheets 7, 8T and 15T + 19T (options *VD-2* and *B-2*) it is possible to achieve 10.5 points, and from the union of datasheets 7, 15T + 19T, and 27T (options *VD-3* and *B-3*) a total of 10 points can be reached.

Therefore, a suitable implementation of the ARERA technical datasheets (that, apart from other things, permits an easy computation of the energy performances of various building and system components) is recommended to be ancillary utilized with the EU criteria in order to assess a unique scheme for the application of the EU Ecolabel brand.

In addition, the above verified correspondence, allowed us to introduce a criterion for ranking the effectiveness of the proposed measures within the framework of the nearly Zero Energy Buildings approach (nearly Zero Energy Agritourism, in this case). In other words, once the ARERA datasheets have provided useful energy saving results, achieved thanks to the implementation of the proposed interventions, it is easy to report such results in terms of closeness to a zero energy situation for a given agritourism.

#### **5. Conclusions**

Agritourisms represent an important reality in the Italian tourism sector, specifically in Sicily due to their numerical consistency and constantly growing trend. The idea at the base of the presented work stems from some considerations regarding the use of a simple method, based on the ARERA technical datasheets (which constitute an official Italian reference), to assess the energy, environmental, and economic benefits related to the implementation of some energy efficiency measures on a given agritourism, specifically aimed at achieving the EU Ecolabel environmental excellence brand, in the perspective of approaching a potential nearly Zero Energy condition.

The results of the conducted analysis put in evidence some discrepancies regarding the application of the ARERA calculation methods, devised for the residential sector, in a wider context, like that of agritourism. Such an outcome was foreseeable, but it has probably been highlighted even more by the fact that the datasheets results are outdated, having not been updated in the last few years.

Nevertheless, the adoption of the proposed efficiency interventions, despite not being specifically defined for the agritourism context, contributed in addressing both structures toward a nearly Zero Energy path, hence, improving their performances in terms of sustainability.

Apart from the interventions proposed by the ARERA, clearly agritourism sites can be interested in further renewable technologies in order to promote their energy sustainability. In fact, the application of solutions like micro wind turbines, biomass, and high efficiency cogeneration for such purposes has been demonstrated [73]. Similar and/or recently available technologies could represent a driver for implementing new ARERA technical datasheets, in order to render them more compliant with the agritourism context, and the EU targets for energy efficiency and emissions reductions in the civil sector.

In conclusion, it arose that, although it was possible to combine the ARERA technical datasheets with the EU Ecolabel criteria, in order to apply the proposed analysis methodology to the agritourism context in a more efficient way, the existing ARERA technical datasheets should be suitably updated and/or replaced by other more effective tools, expressly planned for the accommodation and catering business sector.

**Author Contributions:** All authors contributed equally to the design, experimental analyses and editing of this research. Conceptualization, L.C., M.L.G., G.P., G.R., G.S.; Data curation, L.C., M.L.G., G.P., G.R., G.S.; Methodology, L.C., M.L.G., G.P., G.R., G.S.; Writing—original draft, L.C., M.L.G., G.P., G.R., G.S.; Writing—review and editing, L.C., M.L.G., G.P., G.R., G.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was carried out within the research funds provided by the XXXIII Cycle Doctoral Course in Energy and Information Technologies of the University of Palermo.

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

#### **Appendix A**


#### **Table A1.** ARERA technical datasheets.


#### **Table A1.** *Cont*.

**Table A2.** EU Ecolabel for Tourist Accommodation Services criteria.



#### **Table A2.** *Cont*.


**Table A2.** *Cont*.

#### **Appendix B**

In the following the equations relating to the ARERA calculation procedures for each considered technical datasheet are given, in order to define the parameters reported in the following, Table A3.

Datasheet N◦ 5, "Replacement of simple glazing with double glazing", allows obtaining the gross primary energy savings (*RL*) achievable per individual building:

$$RL = RSL \times S\_{window} \left[ \text{to} / \, \text{year} / \, \text{building} \right] \tag{A1}$$

where:


Datasheet N◦ 7, "Use of photovoltaic systems with an electrical power of less than 20 kW", allows obtaining the achievable specific gross primary energy savings (*RSL*) for each reference physical unit (*UFR*), represented by a photovoltaic system with electrical power <20 kW:

$$\text{RSL} = k\mathcal{W}\_p \times h\_{\text{eq}} \times k\_1 \times 0.22 \cdot 10^{-3} \text{ [toe/year]} \tag{A2}$$

#### where:


Datasheet N◦ 8T, "Installation of solar collectors for the production of domestic hot water", allows obtaining the annual shares of primary net energy savings (*RNC*) for each reference physical unit (*UFR*), represented by the opening surface (m<sup>2</sup> ) of the installed collectors:

$$\text{RN}\_{\text{C}} = \text{RSN} \times \text{LFR} \text{ [toe/year]} \tag{A3}$$

where:


Datasheet N◦ 15T, "Installation of outdoor air electric heat pumps instead of boilers in newly built or renovated residential buildings", allows obtaining the annual shares of primary net energy savings (*RNC*) for each reference physical unit (*UFR*), represented by a standard apartment, which in terms of square meters of heated surface corresponds to about 80–90 m<sup>2</sup> :

$$\text{RN}\_{\text{C}} = a \times \text{RSL} \times \text{LFR} \text{ [toe/year]} \tag{A4}$$

where:


Datasheet N◦ 19T, "Installation of high efficiency outdoor air conditioners with cooling capacity lower than 12 kWf", allows obtaining the annual shares of primary net energy savings (*RNC*) for each reference physical unit (*UFR*), represented by 1 kW cooling capacity of the air conditioning system at nominal conditions (expressed in actual installed cooling capacity):

$$\text{RN}\_{\text{C}} = a \times \text{RSL} \times \text{LFR} \left[ \text{toe}/year \right] \tag{A5}$$

where:


Datasheet N◦ 27T, "Installation of electric heat pump for domestic hot water production in new and existing plants", allows obtaining the annual shares of primary net energy savings (*RNC*) for each reference physical unit (*UFR*), represented by an electric heat pump water heater for the production of domestic hot water (expressed in number of units):

$$\text{RN}\_{\text{C}} = a \times \text{RSL} \times \text{LFR [toe/year]} \tag{A6}$$

where:


Datasheet N◦ 6, "Wall and roof insulation", allows obtaining the gross primary energy savings (*RL*) achievable per insulated surface unit (m<sup>2</sup> ):

$$RL = RSL \times S\_{wall-root} \text{ [toe/year/builliding]} \tag{A7}$$

where:


Datasheet N◦ 20T, "Thermal insulation of walls and roofs for summer cooling in domestic and service sectors", allows obtaining the annual shares of primary net energy savings (*RNC*) achievable per m<sup>2</sup> of insulated surface unit (*UFR*):

$$\text{RN}\_{\mathbb{C}} = a \times \text{RSL} \times \text{LFR } \left[ \text{toe/year} \right] \tag{A8}$$

where:



**Table A3.** Equations (A1) to (A8) input data parameters for the two considered agritourisms.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Barriers of Consumer Behavior for the Development of the Circular Economy: Empirical Evidence from Russia**

**Svetlana Ratner 1,2 , Inna Lazanyuk 1,\* , Svetlana Revinova <sup>1</sup> and Konstantin Gomonov <sup>1</sup>**


**Abstract:** This paper contributes to the literature on sustainable consumption by in-depth analysis of the factors affecting the probability of 57 different practices of proenvironmental behavior (PEBs) in Russia. The set of studied PEBs includes not only popular energy-saving and wastemanagement practices but also more circular patterns of plastic consumption, shopping, and city mobility. To study real and potential barriers to greening consumer behavior models, we conducted a survey of 623 respondents using a questionnaire developed based on a comparative analysis of similar studies conducted in other countries. The processing of the survey results was carried out using nonparametric statistics due to the absence of a normal distribution of the sample for most of the studied characteristics. The results of the study revealed that the main barriers to sustainable consumption in Russia are the lack of appropriate infrastructure as well as a lack of knowledge. Infrastructural barriers in some situations makes sustainable consumer behavior impossible or inconvenient (in this case, preference is given to other types of consumption), or in some cases necessitates spending additional time and money (then sustainable consumer behavior is not completely denied but practiced less often).

**Keywords:** proecological behavior; circular economy; environmental management; survey; nonparametric methods

#### **1. Introduction**

Circular Economy (CE) describes an economic system based on business models for reusing, recycling, and recovering materials in the production and consumption of goods, works, and services. The transition to the circular economy is necessary to maximize each process's performance in the life cycle of goods or services. The concept of CE opens new radical changes in consumption patterns and lifestyles. Implementing the idea of a circular economy requires rethinking the value chain, a conscious consumer approach, and using new business models such as sharing platforms to extend the life of a product [1]. This principle of improving the economy implies that the global production system is designed to have a no waste concept. All useful elements taken from the environment can be reused and waste from one production chain is the starting material for building another one. In contrast to the economy of the linear type, this stimulates the consumer to change items of consumption always to replace them with newer ones. When the economy transitions according to the circular type, it is essential to extend the life cycle of products as long as possible by the following cycles: (1) reuse, (2) remanufacturing, (3) recycling, (4) disposal [2]. At the same time, the more a product is in the first cycles, the cheaper it is to produce in general than when the product goes for disposal immediately after use (the traditional linear model of the production system).

**Citation:** Ratner, S.; Lazanyuk, I.; Revinova, S.; Gomonov, K. Barriers of Consumer Behavior for the Development of the Circular Economy: Empirical Evidence from Russia. *Appl. Sci.* **2021**, *11*, 46. https://dx.doi.org/10.3390/ app11010046

Received: 11 November 2020 Accepted: 21 December 2020 Published: 23 December 2020

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

**Copyright:** © 2020 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/).

The European Development Plan for a Circular Economy is one of the best-developed practical guidelines for the transition in today's economic environment to the new, more sustainable patterns of production and consumption. The European plan considers the formation of sustainable consumption as one of the priority areas of the circular economy [3].

The transition to the concept of a circular economy requires not only the restructuring of production chains but also the reformatting of many logistic, informational, and managerial links, as well as a change in consumer behavior models [4,5]. The lack of technology and infrastructure readiness to support consumer behavior patterns can inhibit the introduction of proecological consumer behavior patterns. The sustainability of the ecological infrastructure ("eco" -infrastructure) and the living environment as socioecological subsystems necessary for a person is achieved by the ability to adapt in a changing world. By "eco" -infrastructure, we mean a complex of natural, natural-anthropogenic and artificial objects and systems that provide conditions for preserving the human environment. The "eco" infrastructure includes elements of traditional infrastructure (natural resources, all mining systems, waste disposal systems, energy, transport, etc.). Many elements of traditional infrastructure must be environmentally stable in order to be included as components of the "eco" infrastructure that a person must carry out in accordance with ethics and environmental regulation. However, modern world literature in the field of sustainable development shows that it is cultural and economic barriers that have the most significant impact on which of the practices of proecological behavior are spread and which are not [6–8].

Consumer demand is a powerful incentive to create new and transformative industries to increase environmental and social responsibility. However, consumers' readiness for such a radical change in traditional patterns should be formed gradually with the transition from simple resource conservation and waste management practices to more complex ones.

Information on the barriers of proenvironmental conduct in a different cultural, economic, and professional context is scattered across the body of literature [9–12]. However, as far as we know, no previous research has investigated the problem of the barriers of proenvironmental behavior in the context of circular economy development in Russia. This country is only in the beginning of its way to the circular economy but has a comparatively high level of environmental awareness in some groups of consumers and strong cultural traditions in sustainable consumption of food and clothes [13].

Most scholars recognize that environmental awareness is the first important step in the adoption of proenvironmental behavior [14,15]. At the same time, some authors argue that there is still a gap between environmental awareness and proenvironmental behavior and rely on behavioral cost concept in explaining the obstacles of sustainable consumption [14,16]. Since there are no consistent government policies in supporting sustainable consumption in Russia, it is possible to assume that behavioral cost of proenvironmental behavior is high in this country. This makes Russia an interesting case to study both from empirical and theoretical points of view.

The purpose of this paper is to assess barriers of proenvironmental behavior associated with the transition towards a more circular daily practice of energy and water consumption, waste management, city mobility, and shopping in Russia. The study contributes to the literature by bringing new empirical evidence of consumers' attitude toward different forms of proenvironmental behavior in Russia and the cultural influences on it.

To achieve the study's objectives, we used a face-to-face questionnaire survey. The survey group mostly included students from the Peoples' Friendship University of Russia (Moscow) and the Kuban State University (Krasnodar), who studies in several training areas, whose educational programs include courses in environmental management, environmental engineering, etc. This choice of respondents is explained by the fact that this category of consumers is the most informed and the most flexible in forming consumer behavior patterns.

#### **2. Literature Review**

Experts realize that the circular economy is characterized by a restorative and closed nature. In the literature, there are three key features inherent in a circular economy: first, enhanced control over natural resources and maintaining a sustainable balance of renewable resources to conserve and maintain natural capital at an inexhaustible level; second, optimization of consumption processes by developing and distributing products, components, and materials that meet the highest level of reuse; third, the identification and prevention of negative externalities of current production activities in order to improve the efficiency of economic and ecological systems [17].

For the functioning of the circular economy model, the Ellen MacArthur Foundation identifies several essential components of the transition and functioning of the circular economy of groups of activities [17]: (a) Regenerate, (b) Share, (c) Optimize, (d) Loop, (e) Exchange.

Despite the relatively recent interest in circular economics, recycled materials are being returned to manufacturing processes at much lower levels. If this system could be improved, loss of value, dependence on volatile product markets, reduced resource productivity, and externalities in the form of environmental pollution could be avoided [18]. Germany's waste management system is one of the most advanced globally; waste disposal is carried out safely for people and the environment. However, only about 14% of the raw materials used in the industry come from processing; the rest still come from raw materials. The introduction of a circular economy is primarily an information problem.

Promoting proecological behavior in people's daily lives is critical in the circular economy's direction and the industrial and commercial sectors' efforts. Therefore, this issue is relatively well studied in modern scientific literature. As evidenced by the results of numerous studies, awareness of the existence of environmental problems concern about the state of the environment is characteristic of almost all social strategy of society, regardless of gender, age, education, type of activity, etc. [19–21]. However, people's environmental activity—their ability to abandon their habitual consumer behavior to take the path of more responsible behavior in everyday practices—is in its infancy [6,22,23]. There are many internal and external barriers to taking real measures to prevent them or reduce the negative impact and consequences if people are well aware of environmental problems [24]. These barriers can be due to various factors such as traditional values, lifestyles, policies, infrastructure, environmental circumstances, various cultures, and countries [5,25].

A significant amount of research addresses selecting factors that determine the choice in favor of proecological behavior. In studies on the role of the individual behavioral aspect in sustainable development, several primary determinants are highlighted that can influence an individual's choice in favor of proecological actions. Thus, the article [26] states that a combination of social and personal moral standards positively impacts the efficient use of electricity, water, other resources, ecological consumerism, and the desire to recycle materials. A person's attitudes and beliefs, reflecting their awareness of problems and a desire to change the situation also stimulate ecofriendly behavior, although in a less significant way. In papers [27,28], it was noted that in addition to these factors, there are several more: the so-called behavioral control, which means the ease with which an individual can put his intentions into action (the presence of monetary funds, time, opportunity); the intrinsic usefulness of environmental good and the willingness to pay for its offer. First of all, the last two parameters affect the desire to purchase environmentally friendly goods instead of conventional ones. Sustainability is based on a balanced relationship of the triple bottom line—people, profit, and planet [29].

The article [30] describes a study of proecological behavior in Canada and the United States. As a result of this study, the authors found a positive relationship between life satisfaction and proecological behavior. Life satisfaction as well as proecological behavior was predicted by more social behavior.

There was high interest in the study, which was conducted by the students of Malaysian universities. The authors of the study [31] carried out a personal assessment of proecologi-

cal behavior. They assessed such qualities of a student's character as openness, benevolence, and conscientiousness. It was found that conscientiousness was the only trait that most influenced the proenvironmental behavior (PEB). The multilevel analysis [32] showed that the relationship between self-transcendence/self-development values–proecological behavior was weaker among societies with higher cultural and socioecological constraints (for example, with lower values of self-expression and economic development). These results clarify when values can promote or inhibit proecological behavior and highlight the need to consider human interaction and context in understanding how personal factors translate into proecological behavior. The article [33] examined the impact of education on proecological behavior in Thailand's educational institutions. The authors concluded that the higher the level of education, the more often the respondents are ready to participate in environmental actions. However, education does not have any impact on willingness to pay environmental taxes. The study [34] estimated that people who received environmental training would demonstrate 4.7 times more voluntary proenvironmental behavior than those who did not receive any training. The article [35] examined personal factors in proecological tourism behavior from a gender perspective. This study is based on a sample of 347 golfers from 16 European countries. The results confirmed the relationship between ecological habits, personal ability, and attitude to the environment. It is shown that only the interaction between personal capabilities and externally-oriented habits influence environmental attitudes.

In studies [36,37], they used the European countries' example. The authors showed that, despite evidence of a positive influence of education on individuals' proecological behavior, clear conclusions are impossible due to an ambiguous causal relationship. Moreover, it may be overlooked factors that force individuals to get more education and care for the environment.

The article [38] examines students' knowledge and behavior in ecological civilization. Based on 13,404 questionnaires from 152 universities and using a polynomial logit regression model, the authors investigated the relationship between Chinese university students' cognition and behavior. This article divides students' ecological civilization knowledge into two parts. One of which is common sense knowledge that comes from long-term public environmental and environmental education in China. The other is the knowledge of the ecological civilization's national strategy, which comes from the recent intense and comprehensive political advertising. This article concludes that university students with the latest knowledge are more willing to implement ecological civilizational behavior. This result illustrates the role of political propaganda in China for ecological civilization to influence its citizens' behavior and inspires other countries to promote environmental policy and promote public environmental protection through the media.

The paper [39] examines the role of stakeholder engagement in establishing and strengthening the sustainability culture in a company transitioning toward a circular economy. The authors also emphasize that the government and universities play an active role in promoting sustainability culture.

In addition to education, energy conservation, for example, is influenced by other aspects. Thus, authors from China [40] considered factors affecting consumer behavior with energy-saving. The authors concluded that urban inhabitants are more economical in terms of energy consumption than rural inhabitants. The authors confirmed a more responsible attitude towards the environment due to improved knowledge in energy conservation and more developed consumption habits. Scientists from Sweden [41] carried out a comprehensive analysis of energy conservation factors. The authors identified the three most significant factors: age of people, type of house, and income—with the decisive factor being the type of housing. People living in private houses are more likely to save than people living in apartment buildings. Socioeconomic factors influence savings on heating more than attitude to the environment.

An interesting study was conducted [42] on procrastination on household energysaving behavior. The author concluded that environmental awareness increases energy

savings but only affects actions that do not require additional costs (lowering the house's temperature during the absence).

One of the most common methods for studying consumer proecological behavior barriers is to conduct surveys [8]. Although this method has its significant drawbacks, which are also actively discussed in the literature (see, for example, works [7,43,44], it remains a priority tool in primary research. The goal is to obtain an overall picture of the prevalence of a particular consumer behavior model and identify the reasons for the nonproliferation of other models.

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

#### *3.1. Procedure and Participants*

To study real and potential barriers to greening consumer behavior models, we conducted a survey of 623 respondents using a questionnaire developed based on a comparative analysis of similar studies conducted in other countries [6,25]. The study was launched in early 2019 on a pilot group of students (N = 100) from Kuban State University. On the stage of pilot research, a face-to-face survey method was used. A personal survey in the case of small- and medium-scale gives an opportunity to study of phenomena in a real-life context [45]. After processing the results of the pilot study and finding interesting consistent patterns, the study was extended to different age groups and regions and continued until the May of 2020.

On the second stage of survey a method of responded-driven sampling was used [46], in which one interviewer from the invited group (students) interviews several (up to 10) people at a time. It allowed us to increase the speed of the study and the size of the sample.

The procedure of responded-driven sampling led us to the following distribution of respondents by age: the majority of respondents belong to the age category from 20 to 29 years old (55%), besides 8% of respondents belongs to the category from 16 to 20 years old. This makes it possible to test the hypothesis about whether the educational process affects the respondent's values and their attitude to the practices of proecological behavior. In addition, the sample has a reasonably well-represented age group from 30 to 39 (11%), from 40 to 49 (11%) and 50 to 59 (9%). Age categories over 60 are less represented (6%).

The distribution of respondents by involvement in the educational process is the following: 62% of respondents are high school or college students, as well as undergraduate and graduate university students. All other respondents (38%) form the category, which was called "not a student" in order to highlight that they are not involved in the education process at the time.

The distribution of respondents by economic activity is the following: 93% of respondents are economically active and only 7% are not. All working respondents are classified as economically active, regardless of whether they are employed (be it self-employed, sole proprietor, an employee at the enterprise) or students. The economically inactive group includes pensioners, housewives, persons on parental leave, or persons with disabilities. This division of respondents into groups was performed to test the hypothesis of how the availability of free time (it is believed that economically active people have less free time) affects the frequency of using various practices of proenvironmental behavior.

The distribution of respondents by region of residence is uneven. The majority live in Moscow (69%). A significant proportion of respondents live in Krasnodar (17%) and the Krasnodar region (10%). Other regions of Russia are slightly represented in the sample, which is a consequence of the research technique—the main groups were formed precisely in Moscow and Krasnodar based on two large universities. Nevertheless, 2% of respondents are representatives of other regions of the Southern Federal District (except for the Krasnodar region), 1% of respondents live in St. Petersburg, and another 1% live in other regions (Samara Oblast, Kamchatka Krai, Yuzhno-Sakhalin Oblast, etc.).

The level of education of the interviewed respondents in general can be defined as high. The majority of respondents (62%) have higher education (31%—bachelor's level and 41% specialist or master's level). Another 6% of respondents have postgraduate educationpostgraduate, doctoral, or residency). Sixteen percent of respondents have secondary education (6%—general and 10%—professional) and another 16% are university students.

Regarding the distribution of respondents by income level, the majority of respondents (60%) determined their income level using the following description "In general ok, but sometimes I have to save". The share of those who defined their financial situation with the following phrase "I can satisfy all my needs and the needs of my family" is also quite large in the sample (23%). Only 1% of the respondents defined their financial situation as "I live in poverty", 11% chose the phrase "I have to save" to describe their financial situation. Five percent of respondents found it difficult to determine their financial situation or refused to answer.

The distribution of respondents by the level of participation in environmental activities is the following: 39% of respondents participated in environmental campaigns to plant trees, clean garbage, collect wastepaper, glass containers, etc. at least once. Almost the same number of respondents did not participate in any environmental activities (38%). Six percent of respondents filed complaints about any environmental pollution, participated in collecting signatures for appeals to the authorities with demands to improve the environmental situation and another 6% made donations for environmental activities. Only 3% of respondents took part in environmental protests and only 2% participated in two or more environmental events.

#### *3.2. Measure of Attitude to Proenvironmental Behavior*

In this study, we use the concept of ecological behavior proposed by Stern [14]. He understands behavior as an interactive product of a personal value system and a set of contextual factors. Accordingly, ecological behavior is "a product of the interaction of personality characteristics (internal factors) and contextual (external) factors" [14]. Internal factors are associated with citizens' characteristics and include environmental knowledge, environmental values, motivation, locus of control, and personal responsibility. External factors can be divided into several categories: political, economic, social, and technological. Political factors are related to the level and speed of adoption of environmental laws and standards for sustainable development. The economic factors include the state of the country's economy as a whole, the volume and structure of environmental investments, the use of natural resources, measures to restore natural resources, and the demand for environmentally friendly products. Social factors influence environmental competence through social cultural values and traditions, the propensity to accept innovations, and change consumer behavior. External technological factors include the level of development of technologies for the treatment of discharges and emissions, disposal and recycling of waste, secondary use of raw materials, environmental research, research in the field of alternative energy, etc.

As noted above, the questionnaire consisted of three parts. The first part of the questionnaire has fixed data internal factors of the respondents.

The second part of the questionnaire aimed to determine the respondents' attitudes on environmental responsibility issues (general environmental self-awareness). We asked the respondents to answer whether they believe that they can improve their city's environment and indicate the types of environmental measures they took part in at least once. The third part of the questionnaire aimed to assess the frequency and reasons for applying (or not applying) practices of proenvironmental behavior in energy conservation, water conservation, waste management, and reducing the use of disposable products and mobility.

In compiling a list of proenvironmental behaviors, we used a variation of the questionnaire proposed by Lee, Kurisu and Hanaki [6]. We chose this study as a prototype for the following reasons: firstly, it contains a large list of possible practices of proenvironmental behavior (includes 58 practices); secondly, it has already been used by other scientists as a basis for international comparisons [25], which allows us to further (when conducting a more extensive study) compare the results obtained for Russia with the results of other countries. It should be noted that from the entire list of 58 practices of proenvironmental behavior, some are almost unknown in modern Russian society. Nevertheless, they were

included in the study to identify possible patterns of behavior transmitted from generation to generation as a way of lean housekeeping. Several practices were excluded from the list due to Russia's inability to use a contradiction to cultural or legal norms. Instead, several new energy-saving practices have been added. In general, our set of PEBs is 90% similar to the set of Lee, Kurisu, and Hanaki [6].

The selected 57 patterns in our study were divided into six groups: (1) patterns in the field of energy conservation; (2) patterns in the field of water conservation; (3) patterns in the field of waste management; (4) patterns that can be arbitrarily called "against plastic" (reducing the use of disposable tableware, packaging, etc.); (5) shopping patterns; (6) urban mobility patterns. This division is explained by the fact that some groups of patterns (for example, in the field of energy efficiency) have been actively promoting at the state level for more than ten years, while others (for example, deciding in favor of purchasing more environmentally friendly goods) have not yet been stimulated. The list of patterns, divided into groups, is presented in Table 1

When answering the question of how often respondents practice each of the 57 models of proenvironmental behavior, one could choose the answer options "never", "rarely", "often", or "always", which, when processed, were translated into a scale from 1 to 4. Besides, respondents were asked to choose the reasons for the application or nonapplication of these practices. Among the possible reasons for applying the practice were "Habit", "Laziness", "Money saving", "Sense of duty", "Fashion" and among the reasons for nonapplication— "Laziness", "Time consuming", "No consideration", "Forget", "Nobody doing", "Costly", "There are no conditions for application", "I did not know that it was so necessary." The choice of each of the reasons was encoded as a Boolean variable (0 or 1).


**Table 1.** Proenvironmental patterns of consumers' behavior.

**Table 1.** *Cont.*


Further, the investigation for answers to research questions was carried out using descriptive and nonparametric statistics. Nonparametric tests were used in case the studied variables are not distributed normally and are measured on an ordinal scale. Research questions and methods of verification are summarized in Table 2.

**Table 2.** Research questions and methods.


**Table 2.** *Cont.*


#### **4. Results**

The survey data showed that most Russians are interested in the environment (73%). More than half (57%) have taken part in nature conservation activity at least once. Every third respondent spoke in the affirmative on personal participation to improve the city's environmental situation. The relatively widespread opinion among respondents about the need to develop initiatives to protect the environment and the degree of responsibility for the environment indicates the lack of necessary changes. From the point of view of the locus of control the respondents' opinions were divided in groups, which are quite close to such a factor as responsibility and determines the measure of a person's moral attitude towards other people and the world around him. Citizens with an internal locus of control who recognize each resident's possibility and ability to influence the current environmental situation were 54% of respondents. On the contrary, 42% of respondents believe that an individual cannot influence the situation with the environment, thereby shifting their responsibility to the relevant social institutions, namely, the federal or local authorities. Despite these differences in answers, Russians are interested in environmental consumption and behavior issues and economic and social challenges. The study's main focus was on determining the statistical relationships between age, involvement in the educational process, economic activity, region of residence, education level, income level, and the degree of participation in environmental activities on the frequency of application of proenvironmental behavior practices.

#### *4.1. Statistical Relationship between the Age and Gender of the Respondent and the Frequency of Application of the Practices of Proecological Behavior*

Since some age categories of respondents were less represented in the sample, we used nonparametric statistical methods for processing research results that allow working with small samples [47,48]. Calculation of nonparametric correlation coefficients Spearman (R) and tau-Kendall (K) revealed the presence of some statistically significant relationships between the age of the respondent and his attitude to various practices of proecological behavior (Table 3). As shown from the nonparametric correlation results, the older generation more often practices more than 50% of the main patterns of proecological behavior. The exception is practice P5—using stairs instead of elevators, for which there is a negative correlation with the respondent's age, which is understandable. The highest values of nonparametric correlation coefficients were obtained in practices: P56—Flame adjustment for cooking, P53—Maintaining air pressure of the tire, P43—Not buying unnecessary products, P51—Avoiding overloading the car, P52—Reducing idling of the car, P50—Doing car checks regularly, P25—Collection and delivery of glass containers to appropriate collection points.

Statistically significant differences in the frequency PEBs between men and women were identified with Mann–Whitney test for practices P4 "Putting hot food into refrigerator after cooling", P50 "Doing car checks regularly", P21 "Avoiding throwing away waste cooking oil", P22 "Following garbage rules", P33 "Throwing away kitchen garbage after it has dried", P35 "Using own bag when going shopping", P37 "Not buying over-packaged products", P47 "Using bicycle or walking" and P49 "Joining the one day without car program". All these PEBs are practiced more often by men than by women.


**Table 3.** Correlations (at the p = 0.05 level) between the respondent's age and the frequency of using the practices of proecological behavior.

#### *4.2. Influence of Involvement in the Educational Process on the Frequency of Application of Practices of Proecological Behavior*

We tested whether education can significantly impact proenvironmental behavior, involving more efficient use of natural resources, recycling materials, and green consumerism. Conducting a series of nonparametric Mann–Whitney tests made it possible to reveal statistically significant differences: in the use of specific practices of proenvironmental behavior between respondents involved in the learning process (schoolchildren, students) and those who completed their education (Table 4).

**Table 4.** Statistically significant (at the *p* < 0.05 level) results of the Mann–Whitney test with a grouping variable involvement in the educational process.



**Table 4.** *Cont.*

Concerning the revealed, statistically significant differences between the groups of respondents receiving education and have completed their education, the second group of respondents demonstrate more conscious proecological behavior in absolutely all aspects.

#### *4.3. The Statistical Relationship between the Grouping Variable—Economic Activity and Frequency of Application Practices Proenvironmental Behavior*

Conducting a series of nonparametric Mann–Whitney tests made it possible to identify statistically significant differences in the use of some proecological behavior practices between economically active and economically inactive respondents (Table 5).



The revealed differences in the frequency of the use of some practices by economically active and economically inactive respondents do not confirm the original author's hypothesis that having more free time in the absence of a permanent job can contribute to the more careful handling of all resources and the use of more labor-intensive waste management practices. On the contrary, economically active respondents demonstrate greater environmental awareness. They are more likely to practice collecting and handing over recycling wastepaper, glass containers, e-waste, drying food waste, and reusable shopping bags. In addition, economically active respondents are more competent in handling household appliances to reduce their energy consumption. At the same time, nonworking respondents more often use those practices associated with saving, first of all, financial resources. Besides, they more often use traditional household practices in cooking, which are aimed not so much at saving resources but simply correspond to the usual established order of things.

#### *4.4. The Impact of the Region of Residence on the Frequency of Use of Practices of Proenvironmental Behavior*

Given the strong uneven distribution of respondents by region of residence, to process the results of statistical testing hypotheses related to regions of residence, nonparametric statistics were also used. In addition, we tried to interpret the results as carefully as possible. The influence of the respondent's region of residence on the frequency of application of proecological behavior practices was investigated using a series of nonparametric Kruskal– Wallis tests. Statistically significant values of the χ 2 -statistic are given in Table 6.

**Table 6.** The results of calculating the Kruskal–Wallis statistics test the hypothesis about the influence of the region of residence on the frequency of using proecological behavior practices.



**Table 6.** *Cont.*

Residents of two regions of Moscow and Krasnodar took part in the survey, and the range of answers was quite comprehensive.

#### *4.5. Effect of Educational Level on the Frequency of Usage Practices Proenvironmental Behavior*

To test the hypothesis about the influence of the respondent's education level on his attitude to various proecological practices, we calculated the nonparametric Spearman (R) and tau-Kendall (K) correlation coefficients. As a result, a weak positive statistically significant result was obtained at the *p* = 0.05 level for some nontrivial practices. It should be noted that Practices P23 (Garbage separation), P25 (Collection and delivery of glass containers to appropriate collection points), P27 (Collection and delivery of used batteries, light bulbs to appropriate collection points), P33 (Throwing away kitchen garbage after it has dried) from the Waste Management category, Practices P28 (Using own cup), P34 (Using receptacle instead of plastic bag) from the category "Rejection of Plastic" and Practices P38 (Buying organic products), P42 (Choosing goods with their CO2 emission in mind (carbon footprint)) from the Purchase category are relatively new for the Russian consumer and/or labor-intensive. Thus, we can talk about the positive influence of the respondent's education level on forming a proecological consumer behavior (Table 7).


**Table 7.** Correlations (at the *p* = 0.05 level) between the respondent's educational level and the frequency of using proenvironmental practices.

#### *4.6. The Impact of Income on the Frequency of Use of Practices of Proenvironmental Behavior*

As a result of testing the hypothesis about the influence of the respondent's income level on the frequency of his use of various practices of proecological behavior by calculating the Spearman and Kendall nonparametric correlations, negative correlations were

revealed between the income level and the use of his own bags for shopping (P35), as well as the refusal to use personal vehicles in urban areas (P47, P48). A weak positive correlation was also found between the respondent's income level and his use of such energy saving practices as cleaning the air conditioner filter (Table 8).

**Table 8.** Correlations (at the *p* = 0.05 level) between the respondent's income level and the frequency of using proenvironmental practices.


Thus, we can say that the level of income has practically no effect on the frequency of using proecological patterns of consumer behavior. The exception is urban mobility practices, where there is a negative impact of income on proecological mobility practices.

*4.7. The Influence of the Level of Involvement in Environmental Activities on the Frequency of Use of Practices of Proenvironmental Behavior*

The influence of the level of involvement in environmental activities (which can be interpreted as the level of environmental awareness) on the frequency of applying proenvironmental behavior was also investigated using a series of nonparametric Kruskal– Wallis tests. Statistically significant values χ 2 -statistics are given in Table 9.

**Table 9.** The results of calculating the Kruskal–Wallis statistics to test the hypothesis about the influence of the respondent's level of involvement in environmental activities on the frequency of using proenvironmental practices.


#### *4.8. Results of the Analysis of Descriptive Statistics on Proenvironmental Behavior Practices*

The calculation of descriptive statistics on the respondents' assessments of the frequency of using the proposed patterns of proecological consumer behavior shows that the most popular practices across all 57 were P1 "Avoiding overloading the refrigerator", P2 "Reducing the opening and closing of the refrigerator door", "P4 "Putting hot food into refrigerator after cooling", P7 "Adjusting the temperature of the air conditioner", P8 "Turning off lights in empty rooms", P10 "Turning off the TV when people are not watching", P11 "Using energy saving mode or turning off when not in use", P28 "Using own cup", P31 "Covering the pan with a lid when cooking or boiling water" er" и P5 P44 "Attempting to fix things before buying a replacement", they have a median score of 2.7 to 3.

0

**P4**

**P8**

**P10**

**P11**

**P2**

**P31**

**P1**

**P44**

**P28**

**P7**

**P52**

**P21**

**P42**

**P27**

**P25**

**P33**

**P49**

**P53**

**P57**

**P55**

0.5

1

1.5

2

2.5

3

0

**P4**

**P8**

**P10**

0.5

1

1.5

2

2.5

3

The least popular practices, with a median score of 1.5 to 2.1, are: P21 "Avoiding throwing away waste cooking oil", P25 "Collection and delivery of glass containers to appropriate collection points", P27 "Collection and delivery of used batteries, light bulbs to appropriate collection points", P33 "Throwing away kitchen garbage after it has dried", P34 "Using receptacle instead of plastic bag", P42 "Choosing goods with their CO<sup>2</sup> emission in mind (carbon footprint)", P49 "Joining the one day without car program", P52 "Reducing idling of the car", P53 "Maintaining the air pressure in the tire", P55 "Using the dishwasher" er" и P5 P57 "Using residual heat when cooking on an electric stove" (Figure 1).

**Figure 1.** Average ratings of 10 most and 10 least popular practices of proecological behavior in Russia.

The most consistent respondents from the regions rated the practices: P1 "Avoiding overloading the refrigerator", P2 "Reducing opening and closing the door of the refrigerator", P3 "Using a lower setting in the refrigerator compartment", P17 "Taking short showers", P36 "Reducing use of disposable products", P41 "Buying ecomark-appliances" (the standard deviation is slightly more than 1). The largest scatter of assessments is observed across practices: P57 "Use of residual heat when cooking on an electric stove", P16 "Turning off the water when washing face or brushing teeth", P51 "Avoiding overloading the car", P53 "Maintaining air pressure of tire". Moscow's respondents did not show consistency in applying practices (standard deviation is on average more than 1.0). The most considerable variation among respondents is observed for practices: P53 "Maintaining air pressure of tire", P55 "Using the dishwasher", P56 "Flame control when cooking on a gas stove", P57 "Using residual heat when cooking on an electric stove" (standard deviation is on average more than 1.5). Among the groups of practices, the most popular were energy saving practices (median 2.8, average 2.6), the least popular practices for waste management (median 1.9, average 2.21) (Figure 2).

One of the central principles of ecological behavior is saving resources in everyday life. Research data showed that a significant proportion of citizens (46%) try to save electricity, water, and gas. Regarding people's behavior, which can be called eco-oriented (saving resources in everyday life, separating garbage, following the norms of waste management), we can conclude that Russians rarely adhere to several pro-eco practices simultaneously, but they are used quite often separately. The most frequently chosen reason for nonapplication of proenvironmental practices of handling household waste was the answer option "There are no necessary conditions for application", which suggests that the overwhelming majority of respondents could make them at least more commonly used with appropriate collection points for glass containers, wastepaper, old clothing, used batteries and light bulbs. However, also popular answers about the reasons for nonapplication were "I see no need", "laziness", "waste of time", "forgetting", which indicates an insufficient level

of environmental self-awareness even in one of the most informed and flexible in household skills group of respondents. The identified differences are most likely explained not so much by the respondent's income level as by other factors that may be indirectly related to the income level: cultural traditions and living conditions. More research is needed to understand better the impact of income on the frequency of using a particular behavior model.

**Figure 2.** Most popular groups of practices.

#### *4.9. Box and Whisker Chart Analysis Results for Proenvironmental Behavior Practices*

The revealed differences were analyzed by constructing box and whisker diagrams, reflecting the median, and quartiles of 25%–75% and the maximum and minimum values of each group's respondents' ratings. Differences due to the discrepancy between only the quartiles of 25–75% and the discrepancy between the median estimates of only groups 4, 5, and 6 (the smallest, cannot be generated randomly) were excluded from the meaningful interpretation. In those cases, when the median values of the most numerous groups of respondents (1, 2, and 3) differed among themselves, we attempted to analyze and substantively explain the statistically significant differences revealed. In Figures 3–16 are box and whisker diagrams with the most noticeable differences in the respondents' assessments between groups 1—Moscow, 2—Krasnodar, 3—Krasnodar region.

**Figure 3.** Box and whisker diagrams for practices P2 (Reducing opening and closing the door of the refrigerator) and P4 (Putting hot food into the refrigerator after cooling) from the category "energy saving".

**Figure 4.** Box and whisker diagrams for practice P6 (Cleaning filter of the air conditioner or cleaner) and P7 (Adjusting the temperature of the air conditioner), the category "energy saving".

**Figure 5.** Box and whisker diagrams for practice P8 (Turning off lights in empty rooms) and P10 (Turning off the TV when people are not watching), the category "energy saving".

**Figure 6.** Box and whisker diagrams for practice P12 (Doing ironing collectively) and P13 (Setting a lower shower temperature), the category "energy saving".

**Figure 7.** Box and whisker diagrams for practice P15 (Using toothbrush cup) and P17 (Taking short showers), category— "Water saving".

**Figure 8.** Box and whisker diagrams for practice P18 (category "Water saving") and P21 (category "Waste management").

**Figure 9.** Box and whisker diagrams for practices P22 (Following garbage rules) and P23 (Garbage separation), the category "Waste management".

**Figure 10.** Box and whisker diagrams for practices P25 (Collection and delivery of glass containers to appropriate collection

**Figure 11.** Box and whisker diagrams by practices P28 (Using own cup), P32 (Composting kitchen garbage), from the categories "No plastic" and "Waste management".

**Figure 12.** Box and whisker diagrams by practices P33 (Throwing away kitchen garbage after it has dried), P34 (Using receptacle instead of plastic bag), from the categories "No plastic" and "Waste management".

**Figure 13.** Box and whisker charts by practices P38 (Buying organic products), P39 (Buying recycled goods), from category "Shopping".

**Figure 14.** Box and whisker charts by practices P41 (Buying ecomark-appliances), P42 (Choosing goods with their CO<sup>2</sup> emission in mind (carbon footprint)), from category "Shopping".

**Figure 15.** Box and whisker diagrams by practices P50 (Doing car checks regularly), P51 (Avoiding overloading the car), category "Mobility".

**Figure 16.** Box and whisker diagrams by practices P52 (Reducing idling of the car), P53 (Maintaining air pressure of the tire), category "Mobility".

Control over opening the refrigerator (P2): in Krasnodar, this practice is used less frequently (median 2.0, which corresponds to the answer "rarely") than in Moscow and small towns of the Krasnodar Territory (median 3.0, which corresponds to the answer "often"). Cooling hot food before placing it in the refrigerator (P4): in Krasnodar and Krasnodar region it is used more often (median 4, which corresponds to the answer "always") than in Moscow (median—3, which corresponds to the answer "often") (Figure 3). Air conditioner filter cleaning (P6) is used much more frequently in Moscow (median 3, which corresponds to the answer "often") than in Krasnodar and Krasnodar region (median 2—"rarely"). Air conditioner temperature control (P7) is more often used in large cities—Moscow and Krasnodar (median 3—"often"), and in Krasnodar region-less often (median 2—"rarely") (Figure 4).

Lighting off in empty rooms (P8) is more commonly used in Krasnodar and Krasnodar region (median 4—"always") than in Moscow (median 3—"often"). Switching off the TV after watching (P10) is most often used in Krasnodar (median 4—"always"), while in Moscow and Krasnodar region—less often (median 3—"often") (Figure 5). Selective ironing (P12) is often practiced in large cities—Moscow and Krasnodar (median 3—"often"). In small towns and rural areas of the Krasnodar region, this practice of energy saving is rarely used (median 2—"rarely). Lowering the water temperature when taking a shower/washing/washing dishes (P13) is more often practiced in Moscow (median 3—"often"), while in Krasnodar and Krasnodar region this practice of energy saving is used less often (median 2—"rarely") (Figure 6).

Using a water cup when brushing teeth (P15) is one of Russia's least popular resourcesaving practices. In Moscow (as well as in the regions that are represented by a small number of respondents), it is still used more often (median 2—"rarely") than in Krasnodar and Krasnodar region (median 1—"never"). Shortening of shower time (P17) as a water saving practice is more often used in Moscow and small towns of the Krasnodar region (median 3—"often") than in Krasnodar (median 2—"rarely") (Figure 7).

Washing dishes in a basin (P18) is more commonly used in Moscow (median 2— "rarely") than in Krasnodar and Krasnodar region (median 1—"never"). A similar situation is observed in the practice of recycling vegetable oil (P21). This situation is rather strange since it is logical to assume that such practices are more accessible to apply in the case of living in your own house and/or running a subsidiary farm, which is much more common in small towns of the Krasnodar region than in such a metropolis as Moscow (Figure 8).

Compliance with waste management standards (P22) is much more common in Moscow (median 3—"rare") than in Krasnodar and Krasnodar region (median 2—"rare"). Judging by the assessments of respondents from regions with a low representation, it can be expected that the same situation of ignoring waste management standards is observed in other regions of Russia. The practice of sorting waste (P23) is more often used in Moscow and small towns of the Krasnodar region (median 2—"rarely") than in Krasnodar (median 1—"never") (Figure 9).

Collection and delivery of wastepaper, glass containers, and electronic waste occurs (P25) more often in Moscow (median 2—"rarely") than in Krasnodar and other settlements of the Krasnodar region (median 1—"never") (Figure 10).

The use of one's tableware instead of disposable ones is equally often practiced in Moscow and Krasnodar (median 3—"often") but even more often in small towns of Krasnodar region (median 4—"always"). Composting of kitchen waste is sometimes practiced in Moscow and small towns of the Krasnodar region (median 2—"rare"), in Krasnodar it does not occur at all (median 1—"never") (Figure 11).

Such a practice of handling household waste as throwing away organic waste only after drying (P33) is sometimes observed in Moscow (median 2—"rarely"), in Krasnodar and other points of the Krasnodar region do not occur at all (median 1—"never"). The use of a container instead of a plastic bag for storing food (P34) is also more often used in Moscow (median 2—"rarely") than in the Krasnodar region and Krasnodar (median 1—"never") (Figure 12).

The purchase of organic products (P38) is more common in Moscow and small towns of the Krasnodar region (median 3—"often") than in Krasnodar (median 2—"rare"). The purchase of processed goods (P39) is sometimes found in large cities—Moscow and Krasnodar but not at all in small towns of the Krasnodar region (median 1—"never") (Figure 13).

Purchases of goods subject to ecolabeling (P41) are more likely to occur in Moscow (median 3—"often") than in Krasnodar and Krasnodar region (median 2—"rare"). The purchase of goods taking into account the carbon footprint (P42) is still found only in Moscow (median 2—"rare"), in the Krasnodar region not found at all (Figure 14).

As for the practices that help reduce the negative impact of the car on the environment, the practice of regular technical inspection and avoidance of overloading the car (P50) is sometimes used in Krasnodar and Krasnodar region (median 2—"rarely"), in Moscow never (median 1) (Figure 15). A similar picture is observed with a decrease in the idling of a car and tire pressure maintenance (P52). (Figure 16). In the Krasnodar region and Krasnodar, they are often used (median 2.5-3), in Moscow—never (median 1).

The respondents' distribution by the level of participation in environmental events held by local communities, regional and federal authorities, and their initiative is shown in Figure 16. Of the most interesting differences between the behavior of respondents, the following can be noted: respondents who show their environmental activity in the form of participation in protest actions are less likely than other groups (including respondents who show their activity in other forms and respondents who do not participate in any what environmental actions) use in practice the patterns P4 (Putting hot food into the refrigerator after cooling), P11 (Using energy-saving mode or turning off when not

in use) and P29 (Avoiding overvolume cooking) from the category "energy saving", P16 (Turning off the water when washing face or brushing teeth) from the category "Water use", P23 (Garbage separation) from the category "Waste management", P28 (Using own cup) from the category "Refusal from plastic", P37 (Not buying over-packaged products) and P38 (Buying organic products) from the category "Shopping", P47 (Using bicycle or walking) and P53 (Maintaining air pressure of the tire) from the category "Mobility". For all other practices, where differences in respondents' behavior from different groups are revealed, this group does not differ from most others. Judging by the identified features of consumer behavior, this group of respondents cannot be classified as environmentally conscious. Their participation in protest actions is most likely motivated not by environmental considerations as by a generally negative attitude towards the authorities' actions.

The most active respondents who took part in more than one type of environmental activities demonstrate a more frequent use of practices P4 (Putting hot food into the refrigerator after cooling) and P56 (Flame adjustment for cooking) from the category "energy saving", P16 (Turning off the water when washing face or brushing teeth) from the category "Water use", P28 (Using own cup) and P37 (Not buying over-packaged products) from in the category "Refusal of plastic", P38 (Buying organic products) from the category "Purchases", P53 (Maintaining air pressure of the tire) from the category "Mobility". At the same time, according to other practices, such as P18 (Washing dishes using jugged water) from the category "Water use", P23 (Garbage separation) from the category "Waste management", P34 (Using receptacle instead of plastic bag) from the category "Avoiding plastic", this group of respondents has the lowest or one of the lowest rates of the frequency of use. Such a spread in the frequency of using various practices of the most active respondents' proecological behavior is most likely evidence that not all of the studied practices are realized as essential and useful. Any striking differences in the behavior of groups of respondents who have never taken part in environmental actions at all, who filed complaints or collected signatures under appeals to the authorities and who participated in tree planting, garbage collection, collection of wastepaper, glass containers, etc., was not identified. This may indicate that the respondents took part in environmental campaigns more out of solidarity, communication, charity, etc., rather than based on clear ideas about the severity of environmental problems. Another interpretation of the results obtained can be offered: it is much easier for even environmentally conscious people to take part in a one-time action than change the usual patterns of their daily consumer behavior.

#### **5. Discussion**

In the extensive literature on various proecological behavior issues, we can highlight several articles that consider various factors and barriers that affect environmental responsibility, which implies a more efficient use of natural resources. Thus, the article [27,49] noted that in addition to environmental awareness, the availability of financial resources, time, and opportunities affect the desire to purchase environmental goods. Studies [50–52] demonstrate that time preferences influence the valuation of PEBs even in very different cultural contexts (US and India), and people who have free time are more engaged in energy-saving proenvironmental behavior. Unlike previous research, our study did not reveal the differences in the frequency of some PEBs between the groups of economically active and inactive respondents, which can support our original hypothesis that having more free time in the absence of a permanent job can contribute to more careful handling of all resources and the use of more labor-intensive waste management practices. The exception is the pattern P8 "turning off lights in empty rooms".

On the contrary, economically active respondents demonstrate greater environmental awareness. In a study [37] on the example of European countries, the author showed evidence of a positive influence of education on the proecological behavior of individuals. The same kind of positive correlation was recently found in [53] on the example of Peru for the PEBs, connected with plastic consumption. Our results also showed that involvement in the educational process has a positive effect on the respondent's values and attitude towards

practices. As for the revealed statistically significant differences between the groups of respondents who are receiving education and have completed their education, the second group of respondents demonstrates more conscious proecological behavior in absolutely all aspects.

Statistically significant differences in the frequency of PEBs for the variable "gender" were found only for nine PEBs, of which three PEBs belong to the mobility group and three PEBs belong to the waste-management group. This result is rather unexpected, since in other countries, for example, in Spain [54], France [55], and China [56], women demonstrate more active proecological behavior.

In the study of authors from China [40], it is shown that other aspects affect energy conservation in addition to education. The authors believe that urban residents are more economical in terms of energy consumption than rural residents. In this, the results of our study are quite consistent with the results of surveys of Chinese scientists conducted in large urban agglomerations and presented in the study [11]. Our results showed that concerning interregional differences, we could say that Moscow is a regional leader in infrastructure development and information support for more "advanced" practices of proenvironmental behavior. The application of these practices requires a certain level of environmental knowledge and specific technical equipment. In small towns and rural areas, traditional patterns of proenvironmental behavior are more often used. However, more detailed conclusions about regional differences can be made only after additional research and expansion of other regions of Russia in the sample.

Despite the fact that the obtained results cannot be directly compared with the parallel results of surveys in other countries presented in the literature, since they were obtained in different macro and micro contexts [56], it can be noted that composting kitchen waste is unpopular in Russia as much as in cities like Tokyo, Bangkok, and Seoul [25]. On the contrary, energy-saving practices are the most popular both in Russia and in all three Asian cities examined with the same methodology. In an article [41,57], the authors identified the most significant factors: age, income, and housing type. They showed that people living in private houses are more likely to save than residents of apartment buildings. The results of the study confirmed the effect of age on the frequency of practice. Namely, the older generation more often practices more than 50% of proecological behavior's main patterns. They are more conscious about their consumption, which can be seen in many practices. However, this behavior is not associated with a decrease in economic activity and a decrease in income. On the contrary, economically active respondents are more likely to practice proenvironmental behavior than economically inactive, especially in cases where proenvironmental behavior is more costly than irresponsible consumption. The main reasons for not applying several consumer behavior patterns are either a lack of understanding of their significance or a complete lack of information about the possibility of such consumer behavior.

The selection of methods for studying barriers to consumer prosustainable behavior confirmed significant shortcomings of the questionnaire used in this study. This method is actively discussed in the literature [7,8,43]. However, when conducting primary research to obtain a general picture of the spread of a particular behavior model and identifying the reasons for the nonproliferation of other models, this method, in our opinion, remains a priority tool. It is possible to identify the main reasons for the weak operation of CE in reality. Central to the effective circular economy design is the observation that there is still a problem with information. For example, potential consumers do not know how many years they can use it when buying a supported product. If the customer can receive such information on time, this may affect his choice. There is also a distorted perception among potential customers that recycled materials are usually considered substandard. In primary markets, external effects (pollution of the environment and air) are not taken into account, which leads to unjustified price advantages of primary materials. It is necessary to increase the literacy of citizens. As the key factors that hinder environmentally sound behavior, the respondents highlight the lack of their initiative and specific environmental measures on the part of the government.

#### **6. Conclusions**

This paper contributes to the literature on sustainable consumption by in-depth analysis of the factors affecting the probability of 57 different PEBs in Russia. The advantages of our research are as follows: the set of studied PEBs includes not only popular energy-saving and waste-management practices but also more circular patterns of plastic consumption, shopping, and city mobility.

The study has revealed that involvement in the educational process (being a student) has a positive effect on the respondent's attitude towards the most energy-saving practices. Besides, the respondents who completed their education demonstrate more conscious proenvironmental behavior in all studied areas: energy and water saving, plastic consumption, food consumption, waste management, and sustainable mobility. At the same time, the probability of PEBs practically not related to the level of income and, therefore, is not the result of intention to save money. Economically active respondents demonstrate greater environmental awareness and are more likely to practice waste collection and recycling and reusable shopping bags. In addition, economically active respondents are more competent in handling household appliances to reduce their energy consumption.

Thus, it is possible to confirm the positive externality of education-higher quality proecological behavior. Given the high cost of recycling modern waste, investment in education can be a preventive measure for improving the environment. In the short term, the effect is hardly possible, but for the regions of Russia, raising awareness and education of the population in the long term can significantly affect the level of environmental awareness.

The main reasons for not applying some sustainable behavior patterns are either a lack of understanding of their significance or a complete lack of information about the possibility of such consumer behavior. Waste management practices are the least popular. We can explain that the unpopularity of waste management practice by the underdevelopment of "eco" infrastructure, which is a barrier to the formation of proecological behavior. As a rule, respondents understand which of the patterns is correct and realize its importance, but the lack of infrastructure at the proper level does not allow using this pattern. The study showed that there are inter-regional differences. We can say that Moscow is a leader in infrastructure development and information support for more "advanced" practices of proenvironmental behavior, the use of which requires a certain level of environmental knowledge and specific technical equipment. In small towns and rural areas, traditional patterns of proecological behavior are more commonly used. However, more detailed conclusions about regional differences can be made only after additional research and expansion of other regions of Russia in the sample. The findings can be used in the educational process and social work with young people.

We realize that the obtained results cannot be directly compared with results of surveys in other countries presented in the literature, since they were obtained in different macro and micro contexts. It should be noted that the factors influencing PEB are still not clearly understood outside a high-income country. However, our major findings do not contradict with conclusions of other authors, arguing that the level of education is one of the most important factors for practicing PEBs.

As for the limitations of our study, we can point out a moderate sample size. This has prevented us from studying the regional differences in behaviors in a more detailed way. Besides the sample size limitation, our study has also neglected the role of cultural values and social norms in stimulating proenvironmental patterns of consumer's behavior.

In our questionnaire design, we have provided questions aimed at identifying the main reasons as to why respondents do not practice PEBs. The answers to these questions should have led us to estimates of behavioral cost. However, in reality, the main reason for not using PEBs was the lack of information about its environmental importance. We will continue our study in trying to redesign the questionnaire, collecting more even samples, and using advanced methods, such as structural modelling, for processing the result of survey.

**Author Contributions:** Conceptualization, S.R. (Svetlana Ratner); Data curation, S.R. (Svetlana Ratner); Formal analysis, I.L.; Funding acquisition, I.L.; Investigation, S.R. (Svetlana Ratner); Methodology, S.R. (Svetlana Ratner); Project administration, I.L.; Resources, K.G., S.R. (Svetlana Revinova); Software, S.R. (Svetlana Ratner); Supervision, S.R. (Svetlana Ratner); Validation, S.R. (Svetlana Ratner), K.G. and S.R. (Svetlana Revinova); Visualization, K.G.; Writing—original draft, I.L.; Writing review and editing, S.R. (Svetlana Ratner), K.G., S.R. (Svetlana Revinova) and I.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** The publication has been prepared with the support of the "RUDN University Program 5-100".

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

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

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

#### **References**


*Article*

## **Which Influencing Factors Could Reduce Ecological Consumption? Evidence from 90 Countries for the Time Period 1996–2015**

#### **Shuai Zhang 1,\*, Dajian Zhu <sup>2</sup> , Jiaping Zhang <sup>2</sup> and Lilian Li <sup>2</sup>**


Received: 18 November 2019; Accepted: 14 January 2020; Published: 18 January 2020

**Abstract:** In the "full world" and Anthropocene, global ecological consumption is beyond natural capital's regenerative and absorptive abilities, and ecological consumption of humanity has to be reduced to have an ecologically sustainable future. To achieve the goal of ecological sustainability, influencing factors that could reduce ecological consumption need to be explored. Based on three panel datasets for the time period 1996–2015, this paper estimates the impacts of urbanization, renewable energy consumption, service industries, and internet usage on ecological consumption for all 90 sample countries, the 42 developed countries, and the 48 developing countries. Education and income are taken as control variables in the panel regressions. As a consumption-side indicator, the ecological footprint is selected to measure ecological consumption. The estimations find that (1) urbanization has negative impacts for all sample countries and the developed countries, and it is insignificant for the developing countries, (2) renewable energy consumption and service industries have negative impacts for all of the three samples, and (3) internet usage has lagged negative impacts for all sample countries, and it is an independent and significant force of reducing ecological consumption in the developing countries rather than the developed countries. It is found that there is a positive linear relationship, an inversed U-shaped relationship, and a U-shaped relationship between ecological consumption and income in all sample countries, the developed countries, and the developing countries, respectively. The estimated results provide guidance for evidence-based policymaking on reducing ecological consumption.

**Keywords:** ecological consumption; influencing factors; panel regressions; ecological footprint

#### **1. Introduction**

In the increasingly "full world" and Anthropocene, and according to the paradigm of strong sustainability, the limiting factor to well-being improvement switched from manmade capital to natural capital [1]. The epochal and fundamental change in the pattern of scarcity warns us that the physical stock of natural capital has to be kept constant, only then enabling an ecologically sustainable future. However, the undisputed fact is that ecological consumption of humanity is beyond natural capital's regenerative and absorptive abilities and that we are living off the "principal" of natural capital [2,3]. Humanity stepped over at least four planetary boundaries, i.e., climate change, rate of biodiversity loss, nitrogen cycle, and change in land use [4,5]. Natural capital is being gradually liquidated and degraded, which could cause declines in provision of ecosystem goods and services and have negative impacts on ecosystem stability and resilience [6].

Therefore, to have an ecologically sustainable future, ecological consumption of humanity has to be reduced until it is kept below natural capital's regenerative and absorptive abilities [7,8]. To achieve the goal of ecological sustainability, influencing factors that could reduce ecological consumption need to be explored in detail, which would provide specific guidance for sound and evidence-based policymaking on reducing ecological consumption. Inspired by previous literature, urbanization, renewable energy consumption, service industries, and internet usage are treated as the latent influencing factors. This paper explores the impacts of the influencing factors on ecological consumption at the global level by selecting representative sample countries covering all of the world. More importantly, this paper explores whether the impacts of the influencing factors on ecological consumption are different or not for countries at different development stages, i.e., developed countries and developing countries. Lessons could be drawn and policy implications could be obtained from the possible estimation differences between developed countries and developing countries.

Because exploring how to reduce ecological consumption from the consumption side is more related to individual lifestyles and consumption habits [9], the ecological footprint (EF) is employed as the proxy of ecological consumption. Regressions based on panel datasets are used to estimate the impacts of the influencing factors on ecological consumption, which could minimize estimation bias caused by omitted explanatory variables. In comparison with time series and cross-section regressions, panel regressions are more capable of controlling econometric problems such as serial correlation and heterogeneity [10]. Two important economic–social factors, education and income, are employed as the control variables in the panel regressions. To some extent, the unobserved and unmentioned influencing factors of ecological consumption could be controlled for by education and income. By incorporating income into the panel regressions, whether the "environmental Kuznets curve" (EKC) hypothesis is valid or not could also be estimated.

The developed countries and the developing countries selected as the samples are those with a population larger than one million and 10 million in 2018, respectively, which to some extent guarantees that their empirical data are relatively reliable [11] (classification of developed countries and developing countries is based on the M49 Standard made by the United Nations, which can be seen from the web page of https://unstats.un.org/unsd/methodology/m49/, accessed on 5 October 2019). Furthermore, the development patterns of developed countries with a population of less than one million and developing countries with a population of less than 10 million are more likely to be unstable and more prone to distortion, which may make the empirical estimations biased and unrepresentative. After sorting out all of the data obtained from public sources, a panel dataset of 90 sample countries for the time period 1996–2015 was available for empirical estimations. The EF data of the 90 sample countries in 2016 were also obtained to estimate the lagged impacts of the influencing factors on ecological consumption, which could reduce estimation bias caused by endogeneity and serve as robust checks for the estimation results. In 2018, the population of the 90 sample countries accounted for 83.85% of the total population, which demonstrates that the sample is quite representative of the whole world and that the empirical findings could provide general guidance on reducing ecological consumption. Among the 90 sample countries, there are 42 developed countries and 48 developing countries.

The remaining of this paper is organized as follows: Section 2 introduces the EF and depicts global and national ecological consumption. Section 3 discusses why the four latent influencing factors are chosen and conducts a literature review. Section 4 presents the regression variables, data sources, and econometric framework. Regression estimations of the impacts of the influencing factors and control variables on ecological consumption for all 90 sample countries, the 42 developed countries, and the 48 developing countries are conducted successively in Section 5. Finally, a discussion and conclusions are presented in Section 6.

#### **2. Ecological Footprint and Levels of Ecological Consumption**

Despite some criticisms of its rationale and methodology, the EF is one of the most popular and inclusive indicators of ecological consumption [12,13]. Some authors argued that the EF is now the most widely used indicator in sustainable development research [14]. The EF measures ecological consumption by calculating the area of biologically productive and mutually exclusive land and water that is required to provide the resources a population demands and to absorb the corresponding wastes in a given year [8]. The EF consists of grazing land footprint (providing animal-based food and other animal products), cropland footprint (providing plant-based food and fiber products), fish product footprint (providing fish-based food products), forest product footprint (providing timber and other forest products), carbon footprint (providing carbon uptake land for absorption of anthropogenic carbon dioxide emissions), and built-up land (representing ecological productivity lost due to occupation of physical space for shelter and other infrastructure) [15].

The EF measures humanity's final demand on a wide range of ecological resources and services from the consumption side (EF*consumption* = EF*production* + EF*imports* − EF*exports*) [16]. Therefore, the land and water to be calculated are not only within national borders but also outside national borders. Because EF values vary greatly with consumption behaviors and habits, it is not difficult to understand global ecological impacts of individual daily lives with the use of the EF [17,18]. The EF is a biophysical rather than monetary accounting approach to measuring ecological consumption. The measuring unit of the EF is global hectares (gha) per capita. A global hectare represents an ecologically productive hectare with global average biological productivity. −

Another prominent advantage of the EF is that it has a counterpart, i.e., the biocapacity (BIO), which measures the theoretical maximum capabilities of ecological systems to meet humanity's demands for ecological consumption. The measuring unit of the BIO is also gha per capita. Comparing national EF to globally available BIO provides a quantitative criterion to assess whether national ecological consumption exceeds globally average ecological capacities. For countries, if their EF values are higher than globally available BIO values, they are countries with an ecological deficit; otherwise, they are countries with an ecological surplus.

Data of the EF and BIO were obtained from the National Footprint Account results (2019 Edition) provided by Global Footprint Network (GFN). Figure 1 depicts temporal trends of the global EF and BIO from 1961 to 2016. The globally available BIO values declined continuously from 3.12 gha to 1.63 gha per capita. The global EF values increased from 2.28 gha per capita in 1961 to 2.87 gha per capita in 1973 and fluctuated between 2.54 gha and 2.87 gha per capita for the time period 1973–2016. Following 1970, the global EF values were larger than the globally available BIO values, which demonstrates that humanity's ecological consumption exceeded the regenerative and absorptive capacities of natural capital and that humanity is living in a state of an ecological deficit.

**Figure 1.** Temporal trends of global ecological footprint (EF) and biocapacity (BIO) (1961–2016). Data source: National Footprint Account results (2019 Edition) from the Global Footprint Network (GFN).

By calculating the ratio of global EF to globally available BIO, how many Earths are needed to support humanity to be ecologically sustainable could be obtained. Figure 2 depicts evolution of the "number of Earths" needed. In 1961, 1970, 1980, 1990, 2000, 2010, and 2016, 0.73, 1.00, 1.19, 1.29, 1.38, 1.67, and 1.69 Earths were needed to support humanity to be ecologically sustainable, respectively. Before 1970, one Earth was sufficient to support humanity to be ecologically sustainable. Since 1970, we needed more than one Earth and, in general, more Earths.

**Figure 2.** Evolution of "number of Earths" needed (1961–2016). Data source: National Footprint Account results (2019 Edition) from the GFN.

As the most influential countries, the G20 countries are used as examples to illustrate national ecological consumption. Table 1 lists the EF values and "number of Earths" of the G20 countries from 1990 to 2016 ("number of Earths" is the ratio of national EF to globally available BIO, which means that, if the level of ecological consumption of one certain country is universal, this is how many Earths would be needed to support humanity to be ecologically sustainable). For the time period 1990–2016, only the EF values of India were lower than the globally available BIO, and India was the only country with an ecological surplus; China and Indonesia transformed from countries with an ecological surplus into countries with an ecological deficit; the EF values of the United States and Canada were extremely large. In 2000, the EF value of the United States was even higher than 10 gha per capita and 5.52 Earths were needed to support the lifestyle of the United States.


**Table 1.** Ecological footprint (EF) values and "number of Earths" of G20 countries for the time period 1990–2016.

Notes: The EF value of Russia of 1990 is not available and the data in the row is of the year 1992.

It is encouraging to find that the EF values of four developed countries (France, Germany, Japan, and the United Kingdom) show downward trends. The EF values of Germany decreased by the largest

extent (29.84%). In 2016, the EF value of Germany was about 40% lower than the EF of the United States. The EF values of seven countries (Argentina, China, India, Indonesia, Turkey, Republic of Korea, and Saudi Arabia) show upward trends. The EF values of Republic of Korea, China, and Saudi Arabia increased by 60.63%, 136.72%, and 193.10%, respectively.

#### **3. Influencing Factors and Control Variables**

Why the four influencing factors and two control variables are selected is discussed in detail. The related and most recent literature is reviewed. In terms of the relationships between the four influencing factors and ecological consumption, four hypotheses are established.

#### *3.1. Urbanization*

As an important economic and social transformation, urbanization is treated as one of foremost influencing factors of ecological consumption. However, the impacts of urbanization on ecological consumption are controversial. On the one hand, because urbanization typically goes hand-in-hand with the industrial process, urbanization would increase ecological consumption due to more consumption of fossil fuels, construction land, cars, electric appliances, and so on. On the other hand, urbanization goes in parallel with a high population density and more ecologically oriented institutions, policies, plans, and technologies, which would permit more efficient use of ecological consumption and, thus, reduce ecological consumption.

By employing the EF as the proxy of ecological consumption and based on panel datasets, Danish and Wang confirmed a positive relationship between urbanization and ecological consumption in 11 emerging countries during 1971–2014 [19], Baloch et al. confirmed a positive relationship in 59 Belt and Road countries for the time period 1990–2016 [20], and Wang and Dong confirmed a positive relationship in 14 sub-Saharan Africa countries for the time period 1990–2014 [21].

By employing the EF as the proxy of ecological consumption, some literature confirmed a negative relationship between urbanization and ecological consumption. Based on time series datasets, Nathaniel et al. confirmed a negative relationship in South Africa for the time period 1965–2014 [22], Ahmed and Wang confirmed a negative relationship in India for the time period 1971–2014 [23], and Dogan et al. confirmed a negative relationship in Nigeria for the time period 1971–2013 [24]. Based on a panel dataset during 1990–2013, Balsalobre-Lorente et al. confirmed a negative relationship in MINT countries (Mexico, Indonesia, Nigeria, and Turkey) [25]. Based on three panel datasets for the time period 1975–2007, Charfeddine and Mrabet found that urbanization would reduce ecological consumption in 15 Middle East and North African (MENA) countries, eight oil-exporting MENA countries and seven non-oil-exporting MENA countries [26].

It is expected that the negative impacts of urbanization on ecological consumption are stronger than the positive impacts and would dominate the relationships between urbanization and ecological consumption. High-density, compact, and modern urban lifestyles are more likely to bring lower levels of ecological consumption. Therefore, the following hypothesis is established:

#### **Hypothesis 1 (H1).** *Urbanization would reduce ecological consumption.*

#### *3.2. Renewable Energy Consumption*

By employing the EF to measure ecological consumption, the estimated relationship between renewable energy consumption and ecological consumption is consistent. Based on two time series datasets, Dogan et al. found a negative relationship between renewable energy consumption and ecological consumption in Nigeria and Turkey for the time period 1971–2013 [24]. Based on panel datasets, Alola et al. confirmed a negative relationship in 16 European Union (EU) countries for the time period 1997–2014 [27], Shujah-ur-Rahman et al. confirmed a negative relationship in 16 Central and Eastern European Countries for the time period 1991–2014 [28], Wang and Dong confirmed a negative relationship in 14 sub-Saharan African countries for the time period 1990–2014 [21], Balsalobre-Lorente et al. confirmed a negative relationship in the MINT countries for the time period 1990–2013 [25], and Olanipekun et al. confirmed a negative relationship in eleven Central and West African countries for the time period 1996–2015 [29].

Renewable energy consumption would reduce greenhouse gas emissions and other pollutants, which constitute major parts of ecological consumption [30,31]. Renewable energy consumption is considered as a key option for reduction in fossil-fuel consumption [32,33]. Whether the EKC hypothesis is valid or not is determined by the significance of renewable energy consumption [10], and that increasing the role of renewable energy consumption is a fundamental strategy in decreasing environmental pressures. In addition, renewable energy consumption may make individuals conscious of ecologically friendly behaviors and lifestyles. Therefore, the following hypothesis is established:

#### **Hypothesis 2 (H2).** *Renewable energy consumption would reduce ecological consumption.*

#### *3.3. Service Industries*

Relative to agricultural and industrial industries, service industries are generally less materialand energy-intensive [29,34]. More importantly, service industries could improve technical efficiencies in using ecological consumption [35]. If service industries account for larger proportions of economic output, national ecological consumption is more likely to be reduced. A paradigm shift from material-intensive and energy-intensive industries to service-centered industries is urgently needed to mitigate negative impacts of ecological overshoot and crisis [36]. Therefore, the following hypothesis is established:

#### **Hypothesis 3 (H3).** *Service industries would reduce ecological consumption.*

#### *3.4. Internet Usage*

Internet changes traditional ways of consumption and production and creates a new economic form, i.e., the internet economy. The internet economy has the potential to reduce ecological consumption in two ways. Firstly, the internet economy is much less material- and energy-intensive than traditional industries; secondly and more importantly, the internet economy could improve the efficiency of every sector of the economy in transforming ecological consumption into economic output [37]. In addition, the internet promotes and facilitates the collaborative economy or the sharing economy, which has the potential to reduce ecological consumption [38,39]. Therefore, the following hypothesis is established:

#### **Hypothesis 4 (H4).** *Internet usage would reduce ecological consumption.*

#### *3.5. Control Variables: Education and Income*

Because education embodies a lot of information on economic and social progress such as scientific and technological progress and human capital accumulation, it is an important and essential control variable. Education may reduce ecological consumption by stimulating ecological awareness and increasing pro-ecological practices. Higher education levels would enable individuals to have more access to various scientific information and knowledge to understand complicated environmental issues and identify causes and consequences of ecological crisis [29]. Furthermore, higher levels of education would increase individual willingness to live an ecologically sustainable life such as installing more renewable energy equipment and participating more in recycling activities [23].

The negative impacts of education may be insignificant because of the attitude–behavior gap [40]. In practice, higher education levels and enough information and comprehension of the ecological crisis may not be transformed into ecologically friendly lifestyles and consumption habits. Moreover, individuals with higher levels of education tend to have high levels of living standards, which often demand higher levels of ecological consumption. Therefore, it is unclear whether education would reduce ecological consumption or not.

Because the EKC hypothesis is quite well known and the EF is a consumption-side proxy of ecological consumption, income is the most important control variable. The EKC hypothesis implies that there is an inversed U-shaped relationship between ecological consumption and income. At low levels of income, levels of ecological consumption and income tend to increase simultaneously. When income reaches a threshold point, levels of ecological consumption would decrease along with further increases in income levels.

Based on a panel dataset of 16 Central and Eastern European Countries during 1991–2014, Shujah-ur-Rahman et al. validated an N-shaped relationship between income and the EF [28]. Based on a panel dataset of 26 EU countries for the time period 1990–2013 and employing the sub-footprints of the EF as the indicators of ecological consumption, Aydin et al. revealed that the EKC hypothesis is not valid [17]. However, based on a panel dataset of 16 EU countries for the time period 1997–2014, Alola et al. found that a 1% increase in real GDP would reduce total EF by 0.81%, which supports the EKC [27]. The above literature shows that the estimations of the EKC hypothesis for developed countries are inconsistent.

By employing the EF to measure ecological consumption, estimations of the EKC hypothesis for developing countries are also inconclusive. Based on time series datasets, Ahmed and Wang confirmed the EKC in India for the time period 1971–2014 [23], and Dogan et al. confirmed the EKC in each of the MINT countries during 1971–2013 [24]. Based on panel datasets, Balsalobre-Lorente et al. validated the EKC in the MINT countries for the time period 1990–2013 [25], but Wang and Dong indicated that economic growth and ecological consumption were positively related in 14 sub-Saharan Africa countries for the time period 1990–2014 [21]. Based on three panel datasets during 1975–2007, Charfeddine and Mrabet showed that the EKC hypothesis was valid in 15 MENA countries and eight oil-exporting MENA countries, and that the relationship between economic growth and ecological consumption was U-shaped in seven non-oil-exporting MENA countries [26].

#### **4. Regression Variables, Data Sources, and Econometric Framework**

To conduct empirical estimations, "urban population (% of total population)" (URB) is employed to measure urbanization, "renewable energy consumption (% of total final energy consumption)" (REN) is employed to measure renewable energy consumption, "services, value added (% of GDP)" (SER) is employed to measure service industries, and "individuals using the internet (% of population)" (INT) is employed to measure internet usage. For the control variables, "mean years of schooling (years)" (MYS) is employed to measure education, and "gross national income per capita, PPP (current international \$)" (GNIPC) is employed to measure income. Abbreviations of all of the regression variables are listed in Table 2.


**Table 2.** List of abbreviations of the regression variables.

Data of URB, REN, SER, INT, and GNIPC were obtained from the World Bank Indicators. Data of MYS were obtained from the Human Development Reports of the United Nations Development Program. All of the data used were obtained from public and reliable data sources, which could guarantee that the regression estimations are repeatable and testable.

After sorting out all of the data of the seven variables and based on the criteria of selecting the sample countries, we finally obtained three panel datasets, i.e., all sample countries (90), the developed countries (42), and the developing countries (48), for the time period 1996–2015. Lists of the developed countries and the developing countries are presented in Tables 3 and 4, respectively. Statistical descriptions of the seven variables for all sample countries, the developed countries, and the developing countries are presented in Tables 5–7, respectively. Furthermore, to estimate the lagged impacts of the independent variables on the EF (lagged by one year) and have as many observations as possible, data of the EF in 2016 for all sample countries were obtained.


**Table 3.** List of 42 developed countries.



**Table 5.** Statistical descriptions of the variables for all sample countries (1996–2015). Obs—observations; Min—minimum; Max—maximum.



**Table 6.** Statistical descriptions of the variables for the developed countries (1996–2015).

**Table 7.** Statistical descriptions of the variables for the developing countries (1996–2015).


Based on panel datasets, the ordinary least square (OLS) was employed to conduct the estimations. We employed Equation (1) to estimate the impacts of the influencing factors and control variables on ecological consumption. Equation (2) is the specification of Equation (1).

$$EF = f(\text{URB}, \text{REN}, \text{SER}, \text{INT}, \text{MYS}, \text{GNIPC}). \tag{1}$$

$$\begin{aligned} \text{LF}\_{i,t} &= \alpha + \beta\_1 \text{IIRB}\_{i,t} + \beta\_2 \text{REN}\_{i,t} + \beta\_3 \text{SER}\_{i,t} + \beta\_4 \text{INT}\_{i,t} + \beta\_5 \text{MYS}\_{i,t} \\ &+ \beta\_6 \text{Ln}(\text{GNIPC})\_{i,t} + \beta\_7 \text{Ln}(\text{GNIPC})^2\_{i,t} + \varepsilon\_{i,t} \end{aligned} \tag{2}$$

*EF* was the dependent variable and *URB*, *REN*, *SER*, *INT*, *MYS*, *Ln*(*GNIPC*), and *Ln*(*GNIPC*) <sup>2</sup> were the independent variables. *Ln*(*GNIPC*) is the natural log form of *GNIPC*. *Ln*(*GNIPC*) 2 is the square of *Ln*(*GNIPC*). Because marginal impacts of income on ecological consumption are supposed to be diminished, *Ln*(*GNIPC*) rather than *GNIPC* is used in Equation (2). Relative to *GNIPC*, *Ln*(*GNIPC*) could minimize the potential estimation bias caused by extreme income values. α represents the intercept term. β1, β2, β3, β4, β5, β6, and β<sup>7</sup> represent the slope coefficients of *URB*, *REN*, *SER*, *INT*, *MYS*, *Ln*(*GNIPC*), and *Ln*(*GNIPC*) 2 , respectively. *i* represents the sample countries (cross-section), which indicates the country-specific effects. *t* denotes the time period (years), which indicates the time series effects. ε*i*,*<sup>t</sup>* is the stochastic error term, which captures the impacts of all unobserved variables on *EF*.

According to the hypotheses in Section 3, β1, β2, β3, and β<sup>4</sup> are expected to be negative. We still could not predict whether β<sup>5</sup> (the coefficient of MYS) is expected to be negative. For ecological consumption and income, there is a monotonically increasing linear relationship if β<sup>6</sup> > 0 and β<sup>7</sup> = 0, there is a monotonically decreasing linear relationship if β<sup>6</sup> < 0 and β<sup>7</sup> = 0, there is an inversed U-shaped relationship if β<sup>6</sup> > 0 and β<sup>7</sup> < 0, which validates the EKC hypothesis, and there is a U-shaped relationship if β<sup>6</sup> < 0 and β<sup>7</sup> > 0, which is contrary to the EKC hypothesis. For the inversed U-shaped or U-shaped relationship, it is easy to be calculated that the turning point values of income are exp (−β6/2β7) (the marginal impacts of *Ln*(*GNIPC*) on *EF* are equal to *dEF dLn*(*GNIPC*) = β<sup>6</sup> + 2β7*Ln*(*GNIPC*). At the turning points, the marginal impacts are zero, and the corresponding values of *Ln*(*GNIPC*) are −β6/2β7. Therefore, the corresponding values of *GNIPC* are exp (−β6/2β7)).

Because values of *Ln*(*GNIPC*) are above zero, an inversed U-shaped or U-shaped relationship between ecological consumption and income cannot be validated if β<sup>6</sup> = 0 and β<sup>7</sup> , 0. Furthermore, it could not be concluded that income is not a significant influencing factor of ecological consumption if β<sup>6</sup> = 0 and β<sup>7</sup> = 0. Under the circumstances, a liner relationship rather than an inversed U-shaped or

U-shaped relationship between ecological consumption and income should be explored and estimated. Therefore, we needed to revise Equation (2) and another estimation equation was proposed. The new estimation equation is as follows:

$$\begin{aligned} \text{LF}\_{i,t} &= \alpha + \beta\_1 \text{URB}\_{i,t} + \beta\_2 \text{REN}\_{i,t} + \beta\_3 \text{SER}\_{i,t} + \beta\_4 \text{INT}\_{i,t} + \beta\_5 \text{MYS}\_{i,t} \\ &+ \beta\_6 \text{Ln}(\text{GNPC})\_{i,t} + \varepsilon\_{i,t} \end{aligned} \tag{3}$$

In order to explore the lagged impacts of the independent variables on *EF* (lagged by one year), the following two equations were estimated:

$$\begin{aligned} \text{EF}\_{i,t} &= \alpha + \beta\_1 \text{IIRB}\_{i,t-1} + \beta\_2 \text{REN}\_{i,t-1} + \beta\_3 \text{SER}\_{i,t-1} + \beta\_4 \text{INT}\_{i,t-1} + \beta\_5 \text{MYS}\_{i,t-1} \\ &+ \beta\_6 \text{Ln}(\text{GNIPC})\_{i,t-1} + \beta\_7 \text{Ln}(\text{GNIPC})^2\_{i,t-1} + \varepsilon\_{i,t-1} \end{aligned} \tag{4}$$

$$\begin{aligned} EF\_{i,t} &= \alpha + \beta\_1 \text{IIRB}\_{i,t-1} + \beta\_2 \text{REN}\_{i,t-1} + \beta\_3 \text{SER}\_{i,t-1} + \beta\_4 \text{INT}\_{i,t-1} + \beta\_5 \text{MYS}\_{i,t-1} \\ &+ \beta\_6 \text{Ln}(\text{GNIPC})\_{i,t-1} + \varepsilon\_{i,t-1} \end{aligned} \tag{5}$$

We could present the regression results by mainly estimating Equations (2) and (4). If an inversed U-shaped or U-shaped relationship between ecological consumption and income could not be validated, Equations (3) and (5) were estimated instead to explore the linear relationship. Based on the three panel datasets, this paper follows the subsequent seven regression procedures:


#### **5. Regression Estimation Results**

Estimation results of the impacts of the influencing factors and control variables on ecological consumption for the three samples are presented. For all sample countries, Models (1)–(6) are presented, and Models (5) and (6) should be used to describe the impacts. For the developed countries and the developing countries, Models (1)–(4) are presented, and Models (3) and (4) should be used to describe the impacts.

#### *5.1. All Sample Countries (90)*

Estimation results of the impacts of the influencing factors and control variables on ecological consumption for all sample countries are presented in Table 8. As can be seen from Models (1) and (2), the estimation results based on the country random effects model and the country fixed effects model were different, especially for the impacts of URB and MYS. Therefore, the Hausman test was conducted to select an appropriate model. The result of the Hausman test (the chi-square statistic was significant at the 1% level) shows that the null hypothesis, i.e., the country random effects model is appropriate, was rejected. Therefore, the country fixed effects model was selected.


**Table 8.** Estimations of impacts of influencing factors and control variables for all sample countries.

Notes: \*, \*\*, and \*\*\* denote significance at the 10%, 5%, and 1% levels, respectively; for Models (1) and (2), probability values (Prob.) are reported in parentheses; for Models (3–6), robust standard errors (RSE) are reported in parentheses.

In Model (3), *Ln*(*GNIPC*) was not statistically significant. In Model (4), neither *Ln*(*GNIPC*) nor *Ln*(*GNIPC*) <sup>2</sup> was significant. According to the arguments in Section 4, an inversed U-shaped or U-shaped relationship between ecological consumption and income could not be statistically validated. Therefore, a linear relationship was explored instead. The scatter plot of EF and GNIPC (Figure 3) further demonstrates that a linear relationship was more appropriate. We ought to interpret the estimated relationships between the independent variables and EF based on the Models (5) and (6).

In Models (5) and (6), the estimated coefficients and significant extents of URB, REN, SER, MYS, and *Ln*(*GNIPC*) showed small differences. URB, REN, and SER had statistically significant and negative impacts on EF. *Ln*(*GNIPC*) had statistically significant and positive impacts on EF. MYS was statistically insignificant. The fact that INT was statistically significant in Model (6) and insignificant in Model (5) demonstrates that INT only had significant lagged impacts on EF. Relative to URB, REN, and SER, the lagged impacts of INT were weaker. The variance inflation factor (VIF) values of MYS and INT were 4.69 and 3.22, respectively, which demonstrates that the estimations of MYS and INT were not likely affected by the problems with multicollinearity.

To sum up, for all sample countries, urbanization, renewable energy consumption, and service industries were the significant influencing factors of reducing ecological consumption; internet usage only had lagged negative impacts on ecological consumption; education had no significant impacts on ecological consumption; ecological consumption and income were positively related. The relationship was linear rather than inversely U-shaped, and the EKC hypothesis was not supported.

**Figure 3.** Scatter plot of EF and gross national income per capita (GNIPC) of all sample countries.

#### *5.2. Developed Countries (42)*

Estimation results of the impacts of the influencing factors and control variables on ecological consumption for the developed countries are presented in Table 9. As can be seen from Models (1) and (2), the estimation results based on the country random effects model and the country fixed effects model were different, especially for the impacts of URB and INT. Therefore, the Hausman test was conducted to select an appropriate model. The result of the Hausman test shows that the country fixed effects model should be selected.


**Table 9.** Estimations of impacts of influencing factors and control variables for the developed countries.

Notes: \*, \*\*, and \*\*\* denote significance at the 10%, 5%, and 1% levels, respectively; for Models (1) and (2), probability values (Prob.) are reported in parentheses; for Models (3) and (4), robust standard errors (RSE) are reported in parentheses.

In Models (3) and (4), the estimated coefficients and significant extents of URB, REN, SER, INT, MYS, Ln(GNIPC), and Ln(GNIPC)<sup>2</sup> showed small differences. URB, REN, and SER had statistically significant and negative impacts on EF. Ln(GNIPC) and Ln(GNIPC)<sup>2</sup> had significantly positive and negative impacts on EF, respectively, which validated an inversed U-shaped relationship between ecological consumption and income. According to the estimated coefficients of Ln(GNIPC) and Ln(GNIPC)<sup>2</sup> in Model (3), the turning point value of GNIPC of the inversed U-shaped relationship was 39,014. The scatter plot of EF and GNIPC (Figure 4) further verifies the turning point. For the sample, there were about 15% of the observations with GNIPC values higher than the turning point. INT and MYS were statistically insignificant. The VIF values of INT and MYS were 3.09 and 2.02, respectively, which demonstrates that the estimations of INT and MYS were not likely affected by the problems with multicollinearity.

**Figure 4.** Scatter plot of EF and GNIPC of the developed countries.

To sum up, for the developed countries, urbanization, renewable energy consumption, and service industries were the significant influencing factors of reducing ecological consumption; internet usage and education had no significant impacts on ecological consumption; the relationship between ecological consumption and income was inversely U-shaped, and the EKC was supported.

#### *5.3. Developing Countries (48)*

Estimation results of the impacts of the influencing factors and control variables on ecological consumption for the developing countries are presented in Table 10. As can be seen from Models (1) and (2), the estimation results based on the country random effects model and the country fixed effects model showed some differences. The Hausman test was conducted to select an appropriate model. The result of the Hausman Test shows that the country fixed effects model was more appropriate and should be selected.

− − − In Models (3) and (4), the estimated coefficients and significant extents of URB, REN, SER, INT, MYS, Ln(GNIPC), and Ln(GNIPC)<sup>2</sup> showed small differences. REN, SER, and INT had statistically significant and negative impacts on EF. Ln(GNIPC) and Ln(GNIPC)<sup>2</sup> had significantly negative and positive impacts on EF, respectively, which validated a U-shaped relationship between ecological consumption and income. According to the estimated coefficients of Ln(GNIPC) and Ln(GNIPC)<sup>2</sup> in Model (3), the turning point value of GNIPC of the U-shaped relationship was 706. The scatter plot of EF and GNIPC (Figure 5) further certifies the U-shaped relationship and the turning point. For the sample, there were about 95% of the observations with GNIPC values higher than the turning point.

− − − −

− − − −

− − − −

− − − −

URB and MYS were statistically insignificant. The VIF values of URB and MYS were 3.24 and 2.82, respectively, which demonstrates that the estimations of URB and MYS were not likely affected by the problems with multicollinearity.


**Table 10.** Estimations of impacts of influencing factors and control variables for the developing countries.

Notes: \*, \*\*, and \*\*\* denote significance at the 10%, 5%, and 1% levels, respectively; for Models (1) and (2), probability values (Prob.) are reported in parentheses; for Models (3) and (4), robust standard errors (RSE) are reported in parentheses.

**Figure 5.** Scatter plot of EF and GNIPC of the developing countries.

To sum up, for the developing countries, renewable energy consumption, service industries, and internet usage were the significant influencing factors of reducing ecological consumption; urbanization and education had no significant impacts on ecological consumption; the relationship between ecological consumption and income was U-shaped, and the EKC hypothesis was not supported.

#### **6. Discussion and Conclusions**

Humanity as a whole is living with an ecological deficit, which is widening in general. Reducing ecological consumption is a basic prerequisite and necessary condition of achieving global ecological sustainability. More and more literature explored the influencing factors of ecological consumption from multiple research perspectives. Based on three panel datasets for the time period 1996–2015, this paper adds to the literature by exploring the impacts of urbanization, renewable energy consumption, service industries, and internet usage on ecological consumption and taking education and income as the control variables.

The main contributions of this paper are as follows: (1) this paper employed the EF as the indicator of ecological consumption. The EF tracks ecological consumption in dimensions of both sources and sinks from the consumption side, which enlarges the discussion on ecological sustainability beyond a certain specific domain and is more illuminating for the demand-side policymaking; (2) by selecting representative sample countries covering all of the world, this paper estimated the impacts of the influencing factors on ecological consumption at the global level, which is helpful of summarizing universal laws of reducing ecological consumption; (3) this paper divided all 90 sample countries into 42 developed countries and 48 developing countries and found the different impacts of the influencing factors for the developed countries and the developing countries, especially the different impacts of urbanization and internet usage.

Urbanization would reduce ecological consumption in all sample countries and the developed countries. However, urbanization was not an independent and significant influencing factor of ecological consumption in the developing countries. Based on a panel dataset for the time period 1980–2015, Adams and Acheampong also demonstrated that urbanization had indeterminate and insignificant impacts on carbon emissions in 46 sub-Saharan Africa countries [30]. As discussed in Section 3.1, the negative impacts do not dominate the potential relationships between urbanization and ecological consumption in the developing countries. To further strengthen the negative impacts and weaken the positive impacts of urbanization on ecological consumption, the urban areas should be more ecologically planed (more vertical and compact rather than horizontal and sprawled), and the urban residents ought to have easier and more convenient access to ecologically efficient technologies and public infrastructure such as consumer durables and mass transit. Unplanned, scattered, and disordered urban sprawl in developing countries must be controlled through joint governance at various levels and by different public departments [41].

For all three samples, renewable energy consumption and service industries would reduce ecological consumption. Developing new, reliable, and affordable green and clean technologies to further promote renewable energy consumption and accelerating economic structure transition (less agricultural and industrial industries and more service industries) are useful and effective ways of reducing ecological consumption. By comparing Tables 9 and 10, it could be found that the negative impacts of both renewable energy consumption and service industries in the developed countries were much stronger than those in the developing countries. The much stronger negative impacts could help explain the inversed U-shaped and U-shaped relationships between ecological consumption and income in the developed countries and the developing countries, respectively [10].

Internet usage was an independent and significant force of reducing ecological consumption in the developing countries. Internet not only brings new information, knowledge, technologies, and goods and services markets, but also more ecologically sustainable lifestyles. However, internet usage was not a significant influencing factor of ecological consumption in the developed countries. By employing the same proxy of internet usage and based on a panel dataset, Salahuddin et al. also found that internet usage could not reduce CO<sup>2</sup> emissions in 31 OECD (Organization for Economic Cooperation and Development) countries for the time period 1991–2012 [42]. Therefore, "internet for ecological sustainability" should be a new guidance principle when developed countries design internet businesses and industries.

Education was not an independent and significant influencing factor of ecological consumption for all three samples. The attitude–behavior gap still exists and prevents ecological awareness from being turned into real and specific ecologically friendly actions. The estimated results remind us that the ecological education we need in the "full world" and Anthropocene does not only have to do with "knowing" but also with "doing". It is hoped that individuals with higher education levels would become the pioneers on living and promoting ecologically sustainable lives.

By employing income as the most important control variable, this paper enriches the long-lasting discussions on the EKC hypothesis and validates three kinds of relationships between ecological consumption and income. For all sample countries, income was a significant factor of increasing ecological consumption, and higher income was not the solution to ecological degradation. For the developing countries, the turning point value of GNIPC of the U-shaped curve was 706, and 95% of the observations were on the right side of the U-shaped curve, which means that, after a very early development stage, ecological consumption began increasing along with further increases in income levels. To our delight, the EKC hypothesis was valid in the developed countries. A positive relationship between ecological consumption and income was significantly reversed when GNIPC values approximately reached 39,014. For the developed countries, because only 15% of the observations were on the right side of the inversed U-shaped curve, more timely and effective measures should be taken to decouple income from ecological consumption and accelerate the arrival of the turning point of the inversed U-shaped curve, which would stimulate developing countries to make corresponding changes and help reduce global ecological consumption by much larger extents.

Much more work needs to be conducted to improve our empirical analysis and further deepen the research context. Because the variables in this paper are defined generally, future research should explore different definitions of the variables, which would give us different perspectives of exploring the influencing factors of ecological consumption. The data quality should also be improved. For example, cultural differences may affect the uniformity of the data in different countries, especially in developed and developing countries. A major drawback of the regression estimations is that we did not find instrumental variables of the influencing factors, especially the instrumental variables of urbanization and internet usage, to conduct sensitivity analyses of the estimation results. Although panel regressions were used to conduct the estimations and lagged impacts were estimated to serve as robust checks for the regression results, the causal relationships between the influencing factors and ecological consumption could not be validated. Additionally, because the regression estimations were on a global level, it was difficult to provide specific measures to reduce ecological consumption of individual countries. Future research can try to combine estimations of global and national levels.

More control variables, for example, political institutional quality and inequality, could be included in the panel regressions in order to refine the estimations. More importantly, the estimated results need further explanation. For example, why the negative impacts of renewable energy consumption and service industries are much stronger in developed countries than in developing countries needs to be explored. The reasons and evidence would provide more specific and scientific policy advice with regard to reducing ecological consumption. Finally, to have an ecologically sustainable future, humanity also has to find solutions to reverse the downward trend of globally available BIO and improve global ecological capacity, which provides more research directions in the field of ecological sustainability.

**Author Contributions:** Conceptualization, S.Z. and D.Z.; methodology, S.Z.; software, L.L.; validation, S.Z. and J.Z.; formal analysis, S.Z.; investigation, S.Z.; resources, D.Z.; data curation, J.Z.; writing—original draft preparation, S.Z.; writing—review and editing, S.Z.; visualization, J.Z.; supervision, D.Z.; project administration, D.Z.; funding acquisition, S.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Shanghai "Chen-Guang Project" (18CG20) and Shanghai Summit Discipline in Design (DA19102).

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

#### **Appendix A**

According to the advice of reviewers, whether the regression results were robust to two different sub-samples (GNIPC values above/below median) was tested. The median value of GNIPC was 9990. The total 1728 observations were divided into 877 observations (GNIPC values above median), which covered 60 countries, and 851 observations (GNIPC values below median), which covered 62 countries. We conducted the tests by employing the same methods and procedures. The regression results of the sub-sample (GNIPC values above median) demonstrated that (1) urbanization, renewable energy consumption, and service industries had significant and negative impacts on ecological consumption, (2) internet usage and education were statistically insignificant, and (3) there was an inversed U-shaped relationship between ecological consumption and income (the EKC hypothesis was supported). As can be seen, the regression results were almost the same as those of the 42 developed countries. The regression results of the sub-sample (GNIPC values below median) demonstrated that (1) renewable energy consumption and service industries had significant and negative impacts on ecological consumption, (2) urbanization, internet usage, and education were statistically insignificant, and (3) there was a U-shaped relationship between ecological consumption and income (the EKC hypothesis was not supported). As can be seen, except for internet usage, the regression results of the variables were almost the same as those of the 48 developing countries. In conclusion, the regression results in this paper were robust to the two sub-samples (GNIPC values above/below median) generally. The details can be seen in Table A1.


**Table A1.** Estimations of impacts of influencing factors and control variables for the two sub-samples (GNIPC values above/below median).

Notes: \*, \*\*, and \*\*\* denote significance at the 10%, 5%, and 1% levels, respectively; for Models (3) and (4), robust standard errors (RSE) are reported in parentheses.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Impact of Virtual Water on Sustainable Development in Gansu Province**

**Weixuan Wang <sup>1</sup> , Jan F. Adamowski <sup>2</sup> , Chunfang Liu 3,4, Yujia Liu <sup>1</sup> , Yongkai Zhang <sup>5</sup> , Xueyan Wang <sup>1</sup> , Haohai Su <sup>1</sup> and Jianjun Cao 1,\***


Received: 14 November 2019; Accepted: 10 January 2020; Published: 13 January 2020

**Abstract:** The concept of virtual water, as a new approach for addressing water shortage and safety issues, can be applied to support sustainable development in water-scarce regions. Using the input-output method, the direct and the complete water use coefficients of industries categorized as primary, secondary, or tertiary, and the spatial flow patterns of the inter-provincial trade in the Gansu province region of China, were explored. The results show that in 2007, 2010, and 2012 the direct and complete water use coefficients of the primary industries were the greatest among the three industry categories, with direct water use coefficients of 1545.58, 882.28, and 762.16, respectively, and complete water use coefficients of 1692.22, 1005.38, and 873.44, respectively; whereas, the direct and complete water use coefficient values of the tertiary industry category were the lowest, with direct water use coefficients of 16.65, 7.74, and 66.89 for 2007, 2010, and 2012, respectively, and complete water use coefficients of 65.46, 66.89, and 72.81 for 2007, 2010, and 2012, respectively. In addition, study results suggest that the volume of virtual water supplied to Gasnu province's local industries has decreased annually, while virtual water exports from the province have increased annually, with the primary industry accounting for 95% of virtual water output. Overall, the virtual water of Gansu province in 2010 showed a net output trend, with a total output of 0.506 billion m<sup>3</sup> , while in 2007 and 2012 it showed a net input trend with a total input of 0.104 and 1.235 billion m<sup>3</sup> , respectively. Beijing, Shanghai, Guangdong, Ningxia and other water-scarce areas were the main input, or import source for Gansu's virtual water; during the years studied, these provinces imported more than 50 million m<sup>3</sup> individually. Based on these results, it is clear that under the current structure, virtual water is mainly exported to the well-developed coastal areas and their adjacent provinces or other water-abundant regions. Therefore, Gansu province should (1) adjust the industrial structure and develop water-saving and high-tech industries; (2) adjust the current trade pattern to reduce virtual water output while increasing its input to achieve balanced economic development and water resource security.

**Keywords:** water resources; virtual water trade; input-output method; Gansu province

#### **1. Introduction**

The concept of virtual water was first proposed by Tony Allan in 1993 [1] to refer to the amount of water used in the production of goods and services [2]. Since 2002, the concept has received extensive attention around the world [3,4]. In China, virtual water has been forwarded as a potential approach to safeguard water resources, especially in water-scarce areas such as the northwest [5]. This interest has prompted many empirical studies on the quantification of virtual water in global crop and livestock products (2352 kg of water in 1 kg of crop, and 6333 kg of water in 1 kg of livestock product, respectively), in Netherlands's and Canada's grain products (1000–2000 kg of water in 1 kg of grain), in Canada's beef products (16,000 kg of water in 1 kg of beef), in American computer chip products (16 kg of water in 1 g of a 32-megabyte computer chip) [6–10], and in Brazil's, Chile's and American paper products (1052.8 kg, 1227.3 kg, 3345.8 kg of water in 1 kg of paper, respectively) [11]. Virtual water studies have also been completed for services, for example, the assessment of virtual water in the Spanish tourism industry (9.7 kg of water per € 1 of tourism income in 2004) [12]. Virtual water trade, which refers to the practice of transporting this hidden water from one country to another, has also been studied. For example, studies have examined the amount and direction of virtual water in agri-food products transported from Italy to China [13], and the factors influencing its characteristics [14]. When employed appropriately, virtual water trade can be used to alleviate water shortages in water-scarce regions and to ensure local water security through the import of water-rich products from water-abundant regions [15,16]. Researchers have investigated the links between virtual water, food security [17], and water footprints [18,19], which is different compared to carbon footprints [20,21]. The former means the cumulative virtual water of all goods and services consumed by one individual or one country [22], while the latter means a measure of the total amount of greenhouse gas emission directly and indirectly produced by an activity or accumulated over the life stages of a product [23,24]. As there is a mutual feedback mechanism between water and carbon footprint, an understanding of water footprint can contribute to reflecting carbon footprint and proposing policies on the utilization of energy resources, and to achieving sustainable development [25–27].

A review of the existing empirical studies on virtual water shows that most are large-scale in scope, and focused on economically developed areas. There is a need for studies at the provincial scale, particularly for the less developed, arid inland provinces of China. Gansu province is located in the central region of northwest China and covers a total area of 425,900 km<sup>2</sup> . It is a typical arid and semi-arid water-deficient area with a dry climate, sparse rainfall and severe water resource shortages. In addition, the unequitable distribution of industrial water use, low utilization efficiency, and over-exploitation of groundwater have resulted in a decrease in river water supply and an increase in desertification, with significant negative consequences for regional social and economic development [28].

The current paper analyzes virtual water content, the primary virtual water export industries and their export destinations, and the spatial patterns of virtual water in Gansu province between 2007 and 2012. The results of this study can be employed to support the sustainable use of water resources and the optimization of industry and trade structures in Gansu province.

#### **2. Statistics and Methods**

#### *2.1. Inter-Regional Input-Output Model*

In 1936, American economist W. Leontief proposed the input-output method to describe the relationships between inputs and outputs in all sectors of an economic system [29,30]. To inform the rational use of water resources, researches have applied the regional input-output method to the analysis of water consumption [31,32]. By compiling the resulting input-output tables and establishing mathematical models, researchers can identify the amount of virtual water in a system, the current direction of its movement [33], and calculate the ultimate water consumption and environmental emissions involved in product production [8,34].

The inter-regional input-output model (IO model) was first proposed by Isard [35] to reflect product trade between regions in a more systematic and comprehensive way compared to the single regional input-output model. The calculation of virtual water trade volumes in a certain area (Figure 1, specific meaning of each formula is presented in Table A1) is based on the assumption that an inter-regional input-output model contains n regions, and that each region has m sectors, so that the mathematical structure of the inter-regional input-output model contains m × n linear equations [36].

**Figure 1.** Flow chart of virtual water calculation [36].

#### *2.2. Data Sources and Processing*

Due to technical limitations and other difficulties associated with data acquisition, compiling the inter-regional input-output tables is demanding in terms of time and resources. These challenges can contribute to delays and irregularities in input-output table publication schedules. In China, for example, the initial data acquisition process for the regional virtual water input-output assessment was initiated in 1987 and has since been repeated at 5-year intervals. Within this process, input-output tables are launched every 3-years, however, these are usually presented 2–6 years after data acquisition for the cycle. The current study utilizes the input-output tables for 2007, 2010 and 2012.

Water consumption data for all industries in Gansu province were derived from the Gansu Water Resources Bulletin in 2007, 2010, and 2012 [37–39]. Forty-two industry sectors were classified into primary, secondary and tertiary industries according to the Industry Classification Regulations [40]. Primary industries include agriculture, forestry, animal husbandry and fishery; secondary industries consist of those in the manufacturing and construction category, while the remaining industries were classified as tertiary (Table A2 presents industry specific classifications). The complete and direct water use coefficients, with the former referring to the water demand of the entire economic system from the perspective of product life cycles, and the latter referring to water used in the production of intermediate products [41,42], were calculated and the volume of virtual water exported from Gansu to other provinces was inferred.

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

#### *3.1. Industrial Virtual Water Consumption in Gansu Province between 2007 and 2012*

#### 3.1.1. Water Consumption Coefficients of Various Industries

Across all three industry categories, the complete water use coefficients are higher than the corresponding direct water use coefficients; this is expected as the concept of complete water use includes both direct and indirect water use (Table 1). In comparison to the secondary and tertiary industries, the complete water use coefficient of the primary industry category is much higher, and the difference between the complete water use coefficient and that of its direct water use coefficient is relatively insignificant, because large quantities of water are directly consumed by many of the industries in the primary industry and their production processes requiring few resources from the secondary and tertiary industries [40]. The secondary and tertiary industries complete water use is significantly higher than direct water use due to the high indirect water use of many of the composite industries that require considerable indirect water use [43,44]. The catering and lighting industries, for example, consume large amounts of agricultural and electrical products, respectively, resulting in large amounts of indirect water consumption during production. Consequently, these industries are labeled the "invisible water-consuming industries" [45,46]. In comparison with the tertiary industry class, the direct water use coefficient of the secondary industry category is relatively large. Secondary industries like the steel and electricity industries consume large amounts of direct water for cooling and rinsing [47,48]. These results suggest that improving water use efficiency in primary industries, and decreasing the use of primary industry products in secondary and tertiary industries, be employed to effectively reduce water consumption across all three industry classes [49]. From 2007 to 2012, the complete and direct water use coefficients of both the primary and secondary industries showed a trend of continuous decline. The decline in the complete water use coefficient in the primary industry is mainly related to the decline in the direct water use coefficient, while for the secondary industry the complete water use coefficient was influenced significantly by both direct water use and the use of water-abundant products of the primary industry [40,50]. The requirements of increasing water use efficiency and water recycling in Gansu province by the 11th Five-Year Plan (Outline of the Eleventh year plan for National Economy and Social Development in the People's Republic of China) from 2006 to 2010 with "six necessaries" principles, including maintaining steady and

rapid economic development, accelerating transformation of economic growth modalities, improving self-directed innovation capabilities, promoting the coordinated development of urban and rural areas, strengthening the construction of a harmonious society, and deepening reforms openings [51], also contributed to the decline of direct water use in the primary and secondary industries. It is estimated that about 43 provincial-level pilot projects were designed to promote water-saving, through optimized and upgraded agricultural and industrial sectors [52]. These measures had significant results. Primarily, Gansu focused on promoting water-saving agricultural techniques specific to local conditions [53]. These techniques contributed to a marked decrease in the agricultural irrigation quota, an increase in grain output, and a significant reduction in the direct water use coefficient of primary industry [54]. The success of promoting appropriate local techniques should be noted - since the impacts of technical and infrastructure changes often vary from region to region, it is suggested that each city in Gansu establish context specific measures towards the promotion of coordinated economic and ecological development province-wide.

**Table 1.** Direct and complete water consumption coefficients of primary, secondary, and tertiary industries in Gansu province from 2007 to 2012.


Unit: m<sup>3</sup> /million yuan.

During the study period, the complete water use coefficient of the tertiary industry class exhibited an increasing trend, while the direct water use coefficient decreased with time. This increase in the complete water use coefficient can be attributed primarily to infrastructure development in the western region of Gasnu Province in response to the 'develop-the-west' strategy under the 11th Five-Year Plan, in which tertiary industry consumed a large number of water-consuming products [48].

3.1.2. Virtual Water Flow and Water Resource Use among Primary, Secondary, and Tertiary Industries

Between 2007 and 2012, the majority of virtual water in Gansu province flowed towards tertiary industry and the smallest portion flowed towards primary industry (Figure 2). Under the 'develop-the-west' strategy of the 11th Five-Year Plan, Gansu focused on the development of infrastructure, basic industries and the tourism belt along the Silk Road. The "One Belt and One Road" has been one of the most important parts of China's strategy of domestic economic and social development, as well as an important part of China's foreign strategy. The "One Belt" refers to the Silk Road Economic Belt, and the "One Road" refers to the 21st Century Maritime Silk Road [55,56]. During this project, the transportation component of the tourism industry has consumed a large amount of water resources to this day [57].

The total volume of water resources demanded by the region is expressed in terms of gross, or total virtual water; this includes locally produced virtual water and that imported from other regions (Table 2) [58]. Compared to 2007, Gansu's primary industry virtual water use decreased by 16.63% in 2010. This reduction is attributed primarily to the emphasis on "grain for green" (conversion of farmland to forests) in the 'develop-the-west' strategy [59], which not only supports reduced virtual water consumption by primary industry, but also affects consumption in the other industry classes that consume raw materials from primary industry. In 2012, the total amount of virtual water increased significantly; this can be ascribed to the continuous development of the national economy which has brought about significant changes in the income level and consumption structure of Gansu residents, and thus contributes to an increase in water resource consumption [48,60]. In addition, the emergence and advancement of the production and service industries has contributed to increased virtual water use across all three industry categories themselves [61].

**Figure 2.** Virtual water flow momentum of primary, secondary, and tertiary industries in Gansu province from 2007 to 2012.

**Table 2.** Total virtual water consumption of primary, secondary, and tertiary industries in Gansu province from 2007 to 2012.


Unit: 100 million m<sup>3</sup> .

#### 3.1.3. Benefit Analysis of Virtual Water in Gansu Province

The results of the current study show a large gap in virtual water use between the three industry categories in Gansu between 2007 and 2012, but the total combined water use did not change significantly during this period. These findings indicate that, although Gansu province has made advances in various industries, inequitable distribution of water resources among the industries is still an issue. Upon examination of both gross production value and virtual water distribution, it appears that water consumption remains stable for each individual industry class, but that gross production value increases, suggesting an improvement in utilization efficiency. The primary industry class consumed the highest proportion of virtual water but exhibited the lowest production value, while the secondary and tertiary industry classes showed the opposite trend (Figure 3).

The added value of primary industries in Gansu province between 2007 and 2012 (17.063, 19.412 and 20.079 billion yuan, respectively) was lower than the added value of the national primary industry class (27.075, 37.895 and 49.391 billion yuan, respectively) [62–67]. This trend is related to grain output and prices; Gansu province's grain output and grain prices were lower than the national level, further indicating that the virtual water distribution in Gansu province in the primary industry is not appropriate. It is suggested that, to remedy this result, Gansu should adjust the distribution of water according to specific conditions (e.g., spatial distribution of water resources within Gansu

province) through trade within the province, encourage the development of water-saving and profitable industries, and reform those industries that are water-consuming and less profitable [68,69].

**Figure 3.** Comparison of the total water production and virtual water consumption in Gansu province from 2007 to 2012.

#### *3.2. Virtual Water Trade in Gansu Province from 2007 to 2012*

#### 3.2.1. Spatial Patterns of Virtual Water Flow in and around Gansu Province

Patterns of virtual water supply and demand between Gansu province and other regions are important to consider. The results of the current study show that, in comparison with 2007, the use of locally produced virtual water in Gansu decreased by 2.262 billion m<sup>3</sup> in 2010 (Table 2). According to the Gansu Provincial Water Resources Bulletin in 2007 and 2010, in 2010 the water-saving irrigated area reached 3 million mu (1/15 ha) more than that in 2007. In 2012, the total amount of virtual water contained in products was 41.259 billion m<sup>3</sup> , and 32.593 billion m<sup>3</sup> of virtual water was provided for local use, accounting for 79.00% of the total. On the basis of 2007 and 2010 statistics, the use of locally produced virtual water increased by 50.72 and 7.334 billion m<sup>3</sup> , respectively in 2012 (Table 2); this is related to the increase in virtual water consumption across all three industry categories. The proportion of virtual water employed locally verses total virtual water decreased annually from 94.8%, 89.52% to 81.43%, while exports have been mounting, from 5.2%, 10.48%, to 18.57% (Table 3). This trend suggests that the increase in total virtual water may be due to a rise in the import of products with high virtual water content, rather than to the decline of virtual water exports. The industries with the highest direct water use are the dominant exporters, such as those in the primary industry category (Table 4). Therefore, Gansu should focus on improving water use efficiency and minimizing the development of industries with high direct water use to alleviate the pressure on water resources.





Between 2007 and 2012, virtual water in Gansu mainly flowed towards developed coastal cities (Figure 4) such as Beijing, Shanghai, Guangdong and Tianjin. This is partially due to the encouragement of advances and reforms in emerging and traditional industries, and the implementation of virtual water strategies in these centers [70,71]. Over 80 million m<sup>3</sup> of virtual water flowed to the provinces of Jiangsu, Zhejiang and Hebei in 2007, 2010, and 2012, respectively. Gansu, as a water-scarce province, consistently exported large volumes of virtual water to water-rich areas from 2007 to 2012. Such a trade pattern no doubt increases the pressure on water resources and affects economic development in Gansu [72,73]. The spatial pattern of virtual water trade in Gansu province was developed by ranking the output volume of virtual water flowing from Gansu to different export regions. The export area included China's eastern coast along with the central, western, and northeast regions. Virtual water output volume was divided into three ranges: less than 80 million m<sup>3</sup> , 40 million to 80 million m<sup>3</sup> , and higher than 40 million m<sup>3</sup> .

**Figure 4.** *Cont*.

**Figure 4.** Spatial pattern of interprovincial flow of virtual water from Gansu province from 2007 to 2012.

3.2.2. Industrial Structure of Virtual Water Inter-Provincial Flow Direction from Gansu Province

The inter-provincial trade of virtual water is carried out using industrial and agricultural products as carriers [74]. From 2007 to 2012, the largest virtual water exporter was primary industry, exporting 1.462 billion m<sup>3</sup> , 2.862 billion m<sup>3</sup> and 7.215 billion m<sup>3</sup> of virtual water in 2007, 2010, and 2012, respectively, accounting for 96.81%, 96.78%, and 97.13% of the total output (Table 4). As such, the virtual water flow of the primary industry determines the trend of the virtual water current in Gansu province as a whole [75]. Because the agricultural sector dominates in Gansu province, it is the main contributor to virtual water exports. From 2007 to 2012, virtual water from primary industries mainly flowed to Beijing, Shanghai, Guangdong, Tianjin, Jiangsu, Zhejiang and Hebei, and the total output volume was over 80 million m<sup>3</sup> , consistently increasing over time. These large export destinations are also the principal virtual water sinks, further supporting the concept that the basic pattern of virtual water output in Gansu province is affected by the trade flow of primary industry products [76]. However, most virtual water produced in Gansu province was consumed there, as agriculture, the primary industry of the province, suffers from low utilization efficiency in terms of irrigation water and rainfall [52,54].

In Gansu province, the secondary industry imported the highest virtual water volume from 2007 to 2012 (Table 4), at 1.294, 1.968 and 5.740 billion m<sup>3</sup> , accounting for 80.17%, 80.26% and 66.24% of the total virtual water import, respectively. This virtual water stream flowed from Jiangsu, Zhejiang, and Henan, each of which retained more than 80 million m<sup>3</sup> during the study time period. This trend increased consistently, as these regions, in the midst of industrialization, consumed large amounts of water for industrial processes [77,78].

3.2.3. Variations and Rationality Analysis of Inter-Provincial Virtual Water Trade in Gansu Province

A net import trend was observed for all virtual water trade in Gansu province in 2007 and 2012; however, in 2010 virtual water displayed a net export trend (Table 4). Gansu province received 104 million m<sup>3</sup> of virtual water in 2007 because the import volume in the secondary and tertiary industry was higher than the export volume of the primary industry. The 'develop-the-west' strategy required the mid-eastern regions to provide facilities, resources, and techniques for the west during that period [55]. The volume of virtual water imported from the primary industry in 2010 was similar

to that of 2007, yet the import volume from the secondary and the tertiary industries increased by 52.09% and 78.47%, respectively.

The export volume of virtual water for the three industry categories combined was two-times greater in 2010 than in 2007; net output of virtual water in 2010 was 506 million m<sup>3</sup> . There are two primary reasons for this; first, although the number of export destinations were reduced, the export volume is extremely large, and second, Gansu increased intermediate inputs and final product exports to other provinces.

Net imports of virtual water were observed in 2012, with a virtual water import volume of 1.235 billion m<sup>3</sup> . The output of the three industry classes combined was twice that of 2010, with the import of the second and tertiary industries, respectively, reaching 2.92 and 7.41 times that of their 2010 counterparts. This is likely related to the policy orientation, industry distribution and state of development of Gansu province, along with the important position of Gansu province in the construction of the "One Belt and One Road" [55].

The distribution of water resources in China is extremely unbalanced. Water resources are more plentiful, and the average annual precipitation is much higher, in coastal areas than in inland areas [79]. The main receivers of Gansu's virtual water exports include Zhejiang, Guangdong and other water-abundant areas, as well as Inner Mongolia, Chongqing, Shaanxi; all regions with more water storage than Gansu [72,80]. Conversely, Ningxia, Beijing, Shanghai and other regions with low water reserve capacity have been exporting water resources to Gansu province. These patterns reflect the regional failure to match the inter-provincial flow of virtual water trade with water source storage.

#### **4. Conclusions and Recommendation**

The current paper systematically explored virtual water flow among industries in Gansu province as well as the inter-provincial virtual water trade patterns. The results of this research indicate that the dynamic variations in virtual water patterns in Gansu province were mediated primarily by the primary and secondary industries from 2007 to 2012. During this period, virtual water for local use decreased, while virtual water exports to developed and water-abundant coastal areas, increased. In addition, it was found that a significant imbalance in virtual water distribution exists across primary, secondary, and tertiary industries, with the majority of virtual water flowing towards tertiary industries, and the smallest volumes flowing towards primary industries.

In view of the imbalance between virtual water distribution across the three industry categories, and the mismatch of inter-provincial virtual water trade with their water resources storage, it is suggested that Gansu emphasizes the promotion of water-saving agricultural techniques, adjusting the structure of the industrial sector by encouraging water-saving and highly-profitable industries, and reducing virtual water export while broadening their import paths, especially from water-abundant areas such as the import of virtual water from adjacent water-abundant areas (Shanxi province and Inner Mongolia), through the establishment of long-term cooperative relationships using appropriate economic policies, to reduce the proportion of water-reliant productions in all industries.

**Author Contributions:** Conceptualization, J.C., W.W. and J.F.A.; Methodology, J.C.; Validation, J.C.; Formal Analysis, W.W.; Investigation, W.W., C.L., Y.L., Y.Z., X.W. and H.S.; Data Curation, W.W.; Writing—Original Draft Preparation, W.W.; Writing—Review and Editing, J.C. and J.F.A.; Visualization, W.W.; Funding Acquisition, J.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Major Program of the Natural Science Foundation of Gansu province, China (18JR4RA002), the Key Laboratory of Ecohydrology of Inland River Basins, Chinese Academy of Science (KLERB-ZS-16-01), the Open Fund for the Key Laboratory of Land Surface Processes and Climate Change in the Cold and Arid Region of the Chinese Academy of Sciences (LPCC2018008), and the Innovation Improvement Project for colleges and universities in Gansu Province in 2019 (2019B-033).

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

#### **Appendix A Calculation of the Direct and Complete Water Use Coe**ffi**cient and Virtual Water Use**


**Table A1.** Calculation of the direct and complete water use coefficient and virtual water use [35,36].

**Table A2.** Industry Specific Classification.


#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Optimal Manufacturing-Reconditioning Decisions in a Reverse Logistic System under Periodic Mandatory Carbon Regulation**

## **Sadok Turki \*, Soulayma Sahraoui, Christophe Sauvey and Nathalie Sauer**

Laboratory of Informatics Engineering and Production of Metz (LGIPM), University of Lorraine, LGIPM, UFR MIM, 3 Rue Augustin Fresnel, F-57070 Metz, France; soulayma.sahraoui9@etu.univ-lorraine.fr (S.S.); christophe.sauvey@univ-lorraine.fr (C.S.); nathalie.sauer@univ-lorraine.fr (N.S.)

**\*** Correspondence: sadok.turki@univ-lorraine.fr

Received: 15 April 2020; Accepted: 18 May 2020; Published: 20 May 2020

**Abstract:** Due to environmental concerns, firms are under increasing pressure to comply with legislations and to take up environmental strategies. This leads researchers and firms to develop new sustainable supply chains, where a new area has emerged for a manufacturing and reconditioning system. The originality of this work consists in simultaneously considering carbon emissions strategies, carbon tax and mandatory emission in a manufacturing-reconditioning system. The proposed system is composed of two parallel machines, a manufacturing stock, a reconditioning stock and a recovery inventory. In order to make the proposed green manufacturing system more realistic, it is assumed that manufactured (new products) and reconditioned products are distinguishable. The quantity of worn products (used products) depends on the sales in the previous periods, and the repair periods of the machines are stochastic and independent. The aim of this work is to determine the optimal capacities of manufacturing and reconditioning stocks that maximize the total profit, as well as the optimal value of worn products under two carbon emissions' limitations. An evolutionary algorithm is developed, along with an efficient improvement method, to find the optimal value of decision variables. Ultimately, numerical results are provided to show the impact of the period of carbon limit and the worn products (returned products) on decision variables.

**Keywords:** production planning; carbon regulation; reconditioning; green logistics; optimization

## **1. Introduction**

Throughout the past few decades, massive carbon emissions have caused serious global environmental damage, such as thick haze and worsening greenhouse gas effects. The Intergovernmental Panel on Climate Change (IPCC), which is the international body for the assessment of climate change, has pointed out in its fifth assessment report that it is necessary to curb the global greenhouse gas (GHG) emissions by 40%–70% from the 2010 level before 2050, and to curb the global GHG emissions to the level of near zero by the end of the 21st century [1]. The IPCC has defined a complete method to standardize the computation of GHG emissions at the national level. Today, most countries are monitoring GHG emissions by using IPCC guidelines to perform annual inventories assessing the quantity of six main gases (CO2: carbon dioxide, CH4: methane, N2O: nitrous oxide, HFCs: hydrofluorocarbons, PFCs: perfluorocarbon, SF6: sulfur hexafluoride) [2]. It was found that CO<sup>2</sup> makes up the broad majority of GHG emissions—it represents about three-quarters of total GHG emissions. For instance, in the USA, the emissions from the industry occupies 24% of total emissions. Indeed, GHG emissions from industry primarily come from burning fossil fuels for energy, as well as GHG emissions from certain chemical reactions necessary to produce goods from raw materials [3,4]. To achieve the IPCC targets,

the development of climate change mitigation technologies have played a pivotal role [5]. Therefore, the industrial companies have adopted new technologies and strategies to alleviate climate change impact associated with their activities. For example, to reduce the carbon emissions, several companies are devoted to recuperating and remanufacturing worn products instead of producing new ones. Indeed, carbon emissions from producing new parts is higher than that for remanufactured ones.

Nowadays, carbon emission control as well as maintaining sustainable economic development has become an increasing challenge, and many countries have attempted to curb carbon emission compulsively through enacting legislation [6]. Several governments promulgated some carbon emission laws, such as carbon taxes, mandatory carbon emissions capacity, and carbon emission cap and trade [7]. The most known regulations are carbon pricing and emissions trading [8]. Those policies can be classified into two categories: mandatory and non-mandatory. The mandatory policy are strict regulations where governments set a strict emissions cap on a firm for a given period. Non-mandatory policies allow companies to choose between reducing and paying for emissions [9], including carbon tax and emission trade (also known as cap and trade). According to the European Commission, European Union Emission Trading System (EU-ETS) is the first and largest emission trading scheme: 75% of international carbon trading covers more than 11,000 companies in 31 countries. Under this policy, companies are allocated a cap or quota on carbon emissions [10]. Indeed, when a company exceeds the allocated cap, it can purchase extra carbon allowance trough the carbon market trading. Conversely, when its carbon emission is under the cap, it can sell its surplus via trading [11]. Besides, carbon tax is a tax on carbon emissions. Companies under this policy are charged a tax proportional to the amount of carbon emitted [12]. Since the implementation of the Kyoto Protocol in 1997, several countries have enacted a variety of carbon tax schemes that have attracted wide attention in the manufacturing management [13]. According to the International Monetary Fund (IMF), carbon tax is the most effective policy to mitigate carbon emissions. With these environmental regulations, companies are constrained to adopt new production policies in order to curb their carbon emissions. From the beginning of the nineties until today, several published papers in the literature have dealt with production decision, taking into account carbon emissions regulations. Ingham and Ulph [14] have presented the case for using a carbon tax to control carbon emissions, and have illustrated the implications for the UK manufacturing. In addition, they have computed a total abatement cost curve corresponding to different levels of emissions' curbing to see whether there are critical points at which costs rise quickly. Fang et al., [15] have developed a mathematical programming model of a flow shop-scheduling problem that considers peak power load, energy consumption, and associated carbon emission in addition to cycle time. The objective is to determine the optimal scheduling that increases the energy consumption and the carbon emissions. Turki and Rezg [16] have proposed an optimal design for a manufacturing/remanufacturing system that sorts worn products into three quality levels. The authors have determined optimal production decisions regarding new and remanufactured products while considering carbon tax policy. As in a real case, the authors have considered that the carbon emissions from producing remanufactured products is different from those of new ones. Indeed, the benefit of the remanufacturing is that the worn products are recovered and reused. In addition, the carbon emission from producing remanufactured products is lower than that for producing new ones. More recently, Dou et al., [17] have considered a manufacturer who produces new parts in the first period and makes new and remanufactured parts in the second period under carbon tax policy, where the tax price differs over the two periods. The authors have determined optimal manufacturing and remanufactured plans that reduce the carbon emissions when the carbon tax varies. He et al., [18] have considered a supply chain network constrained by a stringent mandatory carbon cap, the purpose of which is to examine how stringent carbon regulations and operational decision modes jointly influence the profitability and emissions control of the system. The authors have shown that only when the mandatory cap is set in intermediate level rather than excessively mild or tight can it be effective to balance system profitability and emission control well. Indeed, the mandatory carbon policy is very efficient to curb the carbon emissions; however, due to its sternness, this policy

drives away the investors. Some European countries are periodically applying the mandatory carbon regulation. For example, the French government applies the mandatory carbon policy on transport and manufacturing activities for a defined period when a high pollution rate is revealed. To the best of our knowledge, there is no work in the literature which deals with a periodic mandatory carbon policy. In this paper, we will combine the carbon tax and mandatory carbon regulation. Indeed, we will consider the carbon tax as a permanent policy to reduce the overall emissions, and the mandatory carbon as a periodic regulation when the emissions exceed a tolerated threshold.

One of the processes to curb carbon emissions in the industrial domain is remanufacturing. Remanufacturing, or refurbing, is considered as a great industrial process, within which worn-out products are restored to new products, providing benefits both environmentally and economically, while at the same time enhancing their image as environmentally responsible, since products are reused instead of being discarded, including curbing carbon emissions [19]. It is supposed that remanufactured products are profitable, at 45%–65% of the price of a new product. Furthermore, they are beneficial under emission regulation [20]. In general, the remanufacturing is a series of industrial processes: disassembly of the worn product on parts, restoring of the reusable parts, replacement of the parts if necessary, and reassembly of the parts [21]. Most of the published papers in the literature have assumed that remanufactured products have the same quality as new ones [22–25]. However, in practice, the quality of new products is mostly higher than for remanufactured ones. In this case, the remanufactured parts are called "reconditioned parts" and are sold with a lower price than new ones. Over the lasts years, the reconditioning production keeps increasing. In recent years, for example, in the automotive sector, reconditioning firms cover 35% of all production firms in the sector. That in the aerospace sector represents 25%. Concerning the sales, for example, in Asia, the sales of reconditioned smartphones represent 23% of the mobile phones market. In France and according to BACK MARKET and REMADE TECHNOLOGY companies, which produce and sell remanufactured products such as smart phones, laptops and home appliances, the market of remanufactured products posted a growth rate of 7% in 2018, unlike that of new ones, which fell by 6.5%. In the literature, few researchers distinguish between new and reconditioned parts. Gaur et al., [26] have conducted a real case of a battery manufacturer based in India. The authors have proposed a supply chain system, which considers that reconditioned batteries have a lower performance specification and more limited warranty relative to the equivalent new product. They have determined an optimal sales and production plan as well as configuration of the proposed supply chain system. Turki and Rezg [27] have determined the optimal storage, manufacturing, and reconditioning planning, while taking into consideration the difference between new and reconditioned parts. Indeed, they considered that the reconditioned and new parts are distinguishable, and the reconditioned parts are sold at lower price in a market different from that existing for new ones. Therefore, in order to conduce a real-world case of a manufacturing-reconditioning system, we will consider two distinguished markets: the first market for selling new parts and the second for reconditioned ones. Indeed, we will propose a reverse supply chain system that recovers the worn parts from the first market, and then reconditions them to then be sold in the second market. Furthermore, to satisfy markets demands, we will consider two machines: the first for producing new products and the second for reconditioning worn parts. Moshtagh and Taleizadeh [28] have considered a manufacturing-reconditioning system with two machines, one for manufacturing and one for reconditioning. The authors have taken into account inventory costs, production and reconditioning costs, ordering cost, and sales revenue. However, the authors have neglected some characteristics of manufacturing systems, such as machine breakdowns and time to repair [29]. In fact, in real-life manufacturing systems, the time to repair is a manufacturing dead time that cannot be neglected as it causes production delays and raises the risk of stock shortage. Therefore, when the time to repair is stochastic, as in practice, certainly, the control of the manufacturing, reconditioning, and parts' stocks becomes much more complicated. In this paper, in order to make our proposed system closer to the real-life case, we will consider that both machines are subject to random breakdowns and repairs. Furthermore, to control manufacturing and reconditioning processes, we will apply hedging point policy [30] that ensures that the number of parts does not exceed a defined stock threshold. Moreover, concerning the return of the worn parts, most works have assumed that the return of worn products is proportional to the demands, while this assumption is not true and causes a suboptimal in production policies. Therefore, in this work, we assume that the quantity of worn products depends on the sales in the previous periods. The configuration of the proposed system is inspired by the real-life firm, which produces new and reconditioned products such as smart phones or laptops. In addition, due to already existing environmental preoccupations and hard economic concurrence, we assume that the proposed system is operating in a competitive environment. Thus, in order to keep the system competitive, we aim to find the optimal stocking strategy and production planning that maximize the total profit. This consists to determine the optimal capacities of manufacturing and reconditioning stocks and the optimal production control, under carbon tax and periodic mandatory regulation, taking into account stochastic machine breakdowns and repairs. To perform the study and the simulation of the proposed system closer to reality, we adapt the discrete flow model [31] that is adopted in several works. The benefit of the discrete flow model is that it is faithful to describe discrete production systems. In addition, the simulation of the discrete flow model is clear and not painful. However, its optimization takes a huge amount of time to search the solutions. The efficient optimization method chosen was inspired from the evolutionary algorithm [32], and we developed an improvement method based on local search to give better solutions.

Compared with the existing works, this paper contributes mainly in two ways: the first is to examine a manufacturing-reconditioning system simultaneously considering the periodic mandatory carbon emission and carbon tax. The second is to develop a new approach to determine the quantity of new and reconditioned products to produce when the government imposes a limit carbon emission in a random period, the optimal levels of the manufactured-reconditioned stock, and the optimal percentage of returned worn products that maximizes the total profit. The present work is important for implementation of a new circular economy and green deal, which is now on top of the agenda for Europe.

After having stated our purpose, this paper is organized as follows: Section 2 introduces the manufacturing-reconditioning system and the mathematical models. Section 3 presents the developed optimization method inspired from the evolutionary algorithm to find the optimal value of decision variables. Section 4 analyzes the computational results. Finally, in Section 5, we conclude our work.

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

In this section, we present and explain the proposed system.

In the Table 1 we present parameters and decision variables that are used to formulate the models in this paper:


**Table 1.** Parameters and decision variables.


**Table 1.** *Cont.*

The aim of our work is to find the optimal stock capacity for new products, *SA*\*, for reconditioned products, *SAr*\*, and the percentage of the returned products, *p*\*, that maximizes the total profit, *F*(*t*).

#### *2.1. Description of the System*

Figure 1 presents the model system, that is composed of two machines denoted *M*<sup>1</sup> for manufacturing and *M*<sup>2</sup> for reconditioning, which are subject to random repairs and failures. The machine *M*<sup>1</sup> produces new products from raw material, and we supposed that it never starves. This assumption defines the general case of the manufacturing system where the raw materials are usually acquired. The machine produces new products at variable rate *uA*(*t*) and they are stored in the stock *SA*. The variable *uA*(*t*) represents the decided amount for producing new parts in the period *t*. In fact, *uA*(*t*) is determined at the period *t*-∆*t* according the machine state (θ(*t*)) and the stock level of new parts *sA*(*t*). The production rate *uA*(*t*) takes a value between 0 and its maximum *UA*. The stock *S<sup>A</sup>* satisfies the demand of new products *dA*(*t*) coming from the Market 1. The sold new products will be used for a period α and then will be collected and returned to the manufacturer, then, after, they are stored in the recovery inventory *R*, to be reconditioned by the machine *M*2, and finally are stored in *SAr*. The production rate of the reconditioning machine (*M*2) is denoted *uAr*(*t*), which represents the decided amount for producing reconditioned parts in the period *t*. Indeed, *uAr*(*t*) is determined at the period *t*-∆*t* according the machine state (β(*t*)), the stock level of reconditioned parts *sAr*(*t*) and the stock level of worn products *sAr*(*t*). The production rate *uAr*(*t*) takes a value between 0 and its maximum *UAr*. The stock *SAr* satisfies the demand of reconditioned products *dAr*(*t*) coming from the Market 2. We supposed that the sold reconditioned products in Market 2 will be destroyed after their end of life. Indeed, in the real case, most reconditioned products are reused one time. In order to study the system under the carbon emission regulations used, we take into consideration the carbon emission emitted by producing the new/reconditioned products. In our model, the carbon emissions represent the quantity of carbon emitted from producing either new or reconditioned products. In other words, we do not consider the carbon footprint, but we directly calculate the carbon quantity generated by the production of new and reconditioned parts.

**Figure 1.** Manufacturing-Reconditioning system.

#### *2.2. Mathematical Model*

As assumed before, the machine *M*<sup>1</sup> is always supplied from raw materials, the state of this machine is either up or down. The machine states are presented in the Equation (1).

$$\theta(t) = \begin{cases} 1 & \text{if } \operatorname{machine} M\_1 \text{ is } up \\ 0 & \text{if } \operatorname{ machine} M\_1 \text{ is } down \end{cases} \tag{1}$$

1

The state of machine *M*<sup>2</sup> is given by in the Equation (2).

$$\beta(t) = \begin{cases} 1 & \text{machine } \mathcal{M}\_2 \text{ is } up \\ 0 & \text{machine } \mathcal{M}\_2 \text{ is down} \end{cases} \tag{2}$$

 We supposed that the machines *M*<sup>1</sup> and *M*<sup>2</sup> are subject to a random failure exponentially distributed. When *M*<sup>1</sup> is up, *uA*(*t*) is between 0 and its maximum *U<sup>A</sup>* and when it is down, *uA*(*t*) = 0. It is the same for the machine *M*2. The production rates of the machines *M*<sup>1</sup> and *M*<sup>2</sup> are given by the Equations (3) and (4).

$$\begin{cases} \ 0 < u\_A(t) \le \mathcal{U}\_A & if \theta(t) = 1 \\\ u\_A(t) = 0 & if \theta(t) = 0 \end{cases} \tag{3}$$

$$\begin{cases} \ 0 < u\_{Ar}(t) \le \mathcal{U}\_{Ar} & if \beta(t) = 1 \\\ u\_{Ar}(t) = 0 & if \beta(t) = 0 \end{cases} \tag{4}$$

() 0 () 0 The stocks *sA*(*t*) and *sAr*(*t*) at time *t* are equal to the level of the stock in the previous period plus the number of products produced in the previous period, minus the number of products outgoing from the stock. In order to simplify the calculations, it assumed that ∆*t* = 1. Thus, the stock levels at time *t* are given by the Equations (5) and (6).

$$s\_A(t) = s\_A(t - \Delta t) + u\_A(t - \Delta t) - V\_A(t) \tag{5}$$

*Appl. Sci.* **2020**, *10*, 3534

$$s\_{Ar}(t) = s\_{Ar}(t - \Delta t) + u\_{Ar}(t - \Delta t) - V\_{Ar}(t) \tag{6}$$

The demands *dA*(*t*)/*dAr*(*t*) are satisfied respectively by the stock level in the previous time *sA*(*t* <sup>−</sup> <sup>∆</sup>*t*)/*sAr*(*t* <sup>−</sup> <sup>∆</sup>*t*) if stock levels are higher or equal to *dA*(*t*)/*dAr*(*t*)), otherwise the number of the sold products is equal to the stock levels at *<sup>t</sup>* <sup>−</sup> <sup>∆</sup>*t*. The numbers of new/reconditioned products sold at time *t* are expressed by by the Equations (7) and (8).

$$V\_A(t) = \begin{cases} d\_A(t) & \text{if } s\_A(t - \Delta t) \ge d\_A(t) \\ s\_A(t - \Delta t) & \text{otherwise} \end{cases} \tag{7}$$

$$V\_{Ar}(t) = \begin{cases} d\_A(t) & \text{if } s\_A(t - \Delta t) \ge d\_A(t) \\ s\_A(t - \Delta t) & \text{otherwise} \end{cases} \tag{8}$$

The numbers of unmet demand *PA*(*t*) and *PAr*(*t*) equal the demand at time *t* minus the number of sold products at time *t*. We assume that unmet demands generate a lost cost. This assumption is considered in practice. Indeed, when the customer is not satisfied, a lost cost equivalent to the loss of the business relevance is considered. The Equations (9) and (10) show the number of lost demands at time *t*:

$$P\_A(t) = d\_A(t) - V\_A(t) \tag{9}$$

$$P\_{Ar}(t) = d\_{Ar}(t) - V\_{Ar}(t) \tag{10}$$

The number of worn products that are returned at time *t* is proportional to the quantity of satisfied demand after their lifetime. The number of used products equals a percentage of the number of the products sold at *t* − α, where α is the lifetime of the products and *p* represents the percentage of returned products (0 < *p* < 1). Thus, when *t* − α < 0, there is no return. Therefore, the quantities of returned products are given by the Equation (11).

$$R(t) = \begin{cases} p \cdot V\_A(t - \alpha) & if(t - \alpha) \ge 0 \\ 0 & otherwise \end{cases} \tag{11}$$

The Equation (12) represents the level of the recovery inventory at time *t*, which equals the number of the recovery inventory at *<sup>t</sup>* <sup>−</sup> <sup>∆</sup>*<sup>t</sup>* plus the number of worn products that are returned at time *<sup>t</sup>*, minus the number of products reconditioned by the machine at *<sup>t</sup>* <sup>−</sup> <sup>∆</sup>*t*.

$$s\_r(t) = s\_r(t - \Delta t) + R(t) - \mu\_r(t - \Delta t) \tag{12}$$

For the machine *M*1, the production rate at time *t* (see Equation (13)) depends on the machine state θ(*t*) and on the stock level of new products at time *t*. We assumed that the machine is never starved from raw materials. The production rate follows these constraints:


• Is null when the state of the machine is down.

$$u\_A(t) = \begin{cases} \mathcal{U}\_A & \text{if } \theta(t) = 1 \text{ and } s\_A(t) + \mathcal{U}\_A \le \mathcal{S}\_A\\ \mathcal{S}\_A - s\_A(t) & \text{if } \theta(t) = 1 \text{ and } \mathcal{S}\_A - \mathcal{U}\_A < s\_A(t) < \mathcal{S}\_A\\ d\_A(t) & \text{if } \theta(t) = 1 \text{ or } s\_A(t) = \mathcal{S}\_A\\ 0 & \text{if } \theta(t) = 0 \end{cases} \tag{13}$$

() 1 () ()

For the reconditioned products, the machine *M*<sup>2</sup> is supplied by the recovery inventory. Then, when it is empty or very low, the machine starves. The production rate for reconditioned items looks like the previous one, but also depends on the recovery inventory *R*(*t*). Indeed, the production rate for reconditioned products at time *t* is given the Equation (14).

$$u\_{Ar}(t) = \begin{cases} \mathcal{U}\_{Ar} & \text{if } \beta(t) = 1 \text{ and } s\_{AR}(t) + \mathcal{U}\_{Ar} \le \mathcal{S}\_{Ar} \text{ and } s\_{r}(t) \ge \mathcal{U}\_{Ar} \\ \mathcal{S}\_{Ar} - \mathbf{s}\_{AR}(t) & \text{if } \beta(t) = 1 \text{ and } \mathcal{S}\_{AR} - \mathcal{U}\_{Ar} \le s\_{Ar}(t) < \mathcal{S}\_{AR} \text{ and } s\_{r}(t) \ge \mathcal{S}\_{AR} - s\_{Ar}(t) \\ d\_{Ar}(t) & \text{if } \beta(t) = 1 \text{ and } s\_{AR}(t) = \mathcal{S}\_{Ar} \text{ and } s\_{r}(t) \ge d\_{Ar}(t) \\ \mathcal{S}\_{r}(t) & \text{if } \beta(t) = 1 \text{ and } s\_{AR}(t) + \mathcal{U}\_{Ar} \le \mathcal{S}\_{Ar} \text{ and } s\_{r}(t) < \mathcal{S}\_{Ar} \\ & \text{or } if \beta(t) = 1 \text{ and } \mathcal{S}\_{AR} - \mathcal{U}\_{Ar} \le s\_{Ar}(t) < \mathcal{S}\_{AR} \text{ and } s\_{r}(t) < \mathcal{S}\_{AR} - s\_{Ar}(t) \\ & \text{or } if \beta(t) = 1 \text{ and } s\_{Ar}(t) = \mathcal{S}\_{AR} \text{and } s\_{r}(t) < d\_{AR}(t) \\ 0 & \text{if } \beta(t) = 0 \text{ and } s\_{r}(t) = 0 \end{cases} \tag{14}$$

In this section, we describe the methodology used to plan the production for manufactured and reconditioned products, in the period where the government sets a limitation for carbon emission on a random period (periodic mandatory carbon emission). First, we present the time horizon discretization for the studied problem (see Figure 2).

**Figure 2.** Studied time horizon.

Δ Δ The horizon [0, *T*] is the finite horizon discretized into periods of time ∆*t*, where ∆*t* is the simulation time step. The occurrence of a period under mandatory carbon emission in the system and the length of the period are generated by truncated normal distribution.

• In order to determine the production quantity for manufacturing and reconditioning products under carbon regulation, we have provided a method to determine the quantity that respects the limited carbon emission and answers to the demand according to the determined production quantity under hedging point policy. We recall *qpm* and *qpr* present the carbon emissions of manufacturing and reconditioning respectively, and *q<sup>l</sup>* is the limit quantity of carbon emission to not exceed. We have developed two algorithms (please see Algorithms 1 and 2 that are given in Appendices A and B) that build a function system capable to calculate the proportion of the production rate of the nearest manufactured and reconditioned products under the constraint to not exceed the limit of carbon emissions. The steps are as follows: We calculate the proportion of the production for new and reconditioned products:

$$P\_{\rm SC} = \frac{\mathcal{U}\_A(t)}{\mathcal{U}\_A(t) + \mathcal{U}\_{Ar}(t)} \tag{15}$$

$$P\_{AC} = \frac{\mathcal{U}\_{Ar}(t)}{\mathcal{U}\_A(t) + \mathcal{U}\_{Ar}(t)}\tag{16}$$

• We search, with two counters *i* and *j* running from 1 to the upper bound (*uA*(*t*) and *uAr*(*t*)), to respect the constraint of not exceeding the limit of carbon emission:

$$\text{i.} \, q\_{pm} + \text{j.} \, q\_{pr} < = q\_l \tag{17}$$

• We calculate the proportion:

$$P\_{m2co} = \frac{i}{i+j} \tag{18}$$

$$P\_{R2CO} = \frac{j}{i+j} \tag{19}$$

Then, we developed a function that gives the nearest values corresponding to the sum value *Pm*2*co* and *PR*2*CO* closest to *PSC* and *PAC*. The total profit function is determined by the difference between the total revenue and the total costs. The total revenue includes the total cost from the new products P *T t*=0 *VA*(*t*) · *cvA*, and the total cost for reconditioned products P *T t*=0 *VAr*(*t*) · *cvAr*. The total cost includes: the total production costs for new products P *T T*=0 *uA*(*t*) · *cu<sup>A</sup>* and the total production costs for reconditioned products P *T T*=0 *uAr*(*t*) · *cuAr*, the total storage cost for new products P *T T*=0 *sA*(*t*) · *csA*, the total storage cost for reconditioned products P *T T*=0 *sAr*(*t*) · *csAr*, and the total storage cost for returned products P *T T*=0 *sr*(*t*) · *cs<sup>r</sup>* . The total return cost for worn products P *T T*=0 *R*(*t*) · *cr*, the total lost sales for new products P *T T*=0 *PA*(*t*) · *cpA*, and the total lost sales for reconditioned products P *T T*=0 *PAr*(*t*) · *cpAr*. The total amount of carbon emitted from producing new products over the horizon [0, *T*] is P *T T*=0 *uA*(*t*) · *ct* · *qpm* and from reconditioned products is P *T T*=0 *uAr*(*t*) · *ct* · *qpr*. We try to maximize the total profit function given by:

$$F(t) = \sum\_{t=0}^{t=T} \begin{bmatrix} \left( v\_A(t) \cdot c v\_A + v\_{Ar}(t) \cdot c v\_{Ar} \right) - \left( s\_A(t) \cdot c s\_A + s\_{Ar}(t) \cdot c s\_{Ar} + \right) \\\ s\_r(t) \cdot c s\_r + P\_A(t) \cdot c p\_A + P\_{Ar}(t) \cdot c p\_{Ar} + \mathcal{R}(t) \cdot c r + u\_A(t) \cdot c u\_A + \\\ u\_{Ar}(t) \cdot c u\_{Ar} + u\_{Ar}(t) \cdot c t \cdot q\_{pr} + u\_A(t) \cdot c t \cdot q\_{pm} \end{bmatrix} \tag{20}$$

#### **3. Optimization Method**

The mathematical model presented above is clearly based on simulation. Machines states are simulated, random failures are also simulated according to an exponential distribution law, and time flow is also simulated. From all of this, values of stock levels, number of new and reconditioned products, lost demands, return products quantities, production rates of both machines, limit of carbon emission constraint, production plan, and total profit function are calculated, according to Equations (5) to (20). It is obviously not an integer linear mathematical model, but we have developed an optimization method to cope with the non-linearity and random variables challenges, and managed to obtain a total profit function maximization.

The total profit function was programmed in the language C++ with the free software Dev-C++. The optimization method was also developed in this language. In order to facilitate the method handling, we have put the decision variables in a vector. Each value of the vector can be between a minimum and a maximum value. The optimization process is illustrated in Figure 3.

**Figure 3.** Optimization process of the decision variables (*SA*\*, *SAr*\*, *p*\*)

The model is stochastic by construction. Therefore, we calculate the average of three simulations to obtain a more precise value of the objective function, in order to be able to have the best surface possible to perform optimization. The simulation time horizon is fixed to 10<sup>7</sup> time units. The optimization method is inspired from evolutionary algorithms and returns the optimal values of the vector which maximizes the total profit function. The quality of the results obtained by this method depends on the number of tests. The larger the number of tests, the better the results. However, the simulation time can also increase too much. Thus, we developed an improvement method based on local search that can give better solutions by testing the function the least amount possible. The steps are the following:

Initialization: random individuals are created, with each vector of values respecting the boundaries of each of its parameters.

Evaluation: In addition to the pre-determined constraints on variables, the following constraint is written: *S<sup>A</sup>* > *SAr*. This equation imposes that the storage of new products should be higher than that of reconditioned ones in order to ensure the return of worn products. Consequently, unadapted individuals are rejected.

Selection: We retain the best adapted individuals, those who give satisfying total profit function results.

Mutation: Among the best adapted individuals, we try some neighborhood tests, thanks to a deviation of some values of the vector describing an individual, to see if its characteristics are better adapted.

Evaluation: If a new generated individual meets the constraints, it is evaluated with its total profit value.

Replacement: If the new total profit function is better than the previous, it is retained.

In the end, the final solution is represented by the following vector (*SA*, *SAr*, *p*).

The second part of the proposed optimization method is composed of local improvement procedures, which are very useful to perform a better optimization [32]. We developed three main procedures. The first one changes the value of a parameter, chosen at random, in one sense or in the other, here again chosen at random (positive or negative). It can go up until the appropriate bound, but the variation is also at random (it may not "hit" the bound). The set of parameters is then inside each variable bounds, by construction, but the constraints may be violated anyway. Then, after a constraint checking, its objective function is evaluated. If it improves the current solution, it is kept, and if not, the counter of tested solutions is incremented. We allow to test up to 100 valid neighborhoods.

The second neighborhood iteratively repeats the procedure described below, until there is no improvement between two successive iterations. Positive and negative alterations are tested on each of the chromosomes of an individual. If any best improves the objective function, it is kept. Since this neighborhood is iterative, it is only launched once.

The third neighborhood we have implemented is also iterative. For each couple of values of the decision variables vector, we test an alteration of each possible direction's combinations. Since the vector we have to optimize is composed of 3 values (*SA*\*, *SAr*\*, *p*\*), at each time this neighborhood method is called, it evaluates (*c* 2 3 × 4 = 12) combinations of values around a current vector. A Pseudo-code algorithm of this iterative neighborhood is given in Appendix C (Algorithm 3). The optimization method we have developed combines both the necessity to consider a stochastic function, with the uncertainty and weakness associated, and the need of efficiency with a relatively long computing, time-consuming objective function. As well as the objective function, this method is also stochastic, and its results cannot be considered as optimal. Nevertheless, they allow a quick convergence towards the optimal region by testing an infinitesimal part of the solution space, thanks to the power of the neighborhoods developed. Moreover, in a stochastic solution space, all but convex, and where quasi-optimal solutions have values in a relatively narrow range, the method we have proposed, designed, and developed, optimizes correctly. This allows us first to work on this exhaustive model and to exploit its richness to come to interesting scientific results. This also proves the relevance of this method to treat optimization problems with high computing, time-consuming objective functions.

#### **4. Numerical Results**

This section conducts experimental computations to illustrate the proposed model and explore different results to investigate the impact of interesting parameters on the optimal values of the decision variables. In the beginning of this section, we introduce numerical data, then we present four subsections. In the first, we investigate the impact of the returned products percentage, *p*, on the optimal stock capacity for new products, *SA*\*, and reconditioned products, *SAr*\*. This study analyzes the impact of the quantity of returned worn products on the manufacturing, reconditioning, and stocking decisions. In the second, we examine the influence of the period of mandatory carbon regulation (W) on the optimal stock capacity for new products, *SA*\*, reconditioned products, *SAr*\*, and the percentage of the returned products, *p*\*. This study allows a firm leader to determine the optimal production and stocking planning when the mandatory carbon period changes. In the third, we analyze the influence of the limited quantity of carbon emission imposed by the government on optimal values *SA*\*, *SAr*\*, and *p*\*. This study proposes the optimal stocking management when the government imposes a quantity limit of carbon emissions. In the fourth, we analyze the impact of the carbon emission cost imposed by the government on optimal values *SA*\*, *SAr*\*, and *p*\*.

The used input data are:


The demand *dA*(*t*) is generated by truncated normal distribution, in which the average = 15 and the standard deviation = 7. For *dAr*(*t*), the average = 7 and the standard deviation = 4. The generation of time to repair and time between failures are exponentially distributed. The occurrence of the period of carbon emission limit and length of the period denoted (*W*) are generated by truncated normal distribution. For length of the period *W*, the lower truncated is 40 and the upper truncated is 60, with the average = 50 and standard deviation = 10, and for the occurrence of the period, the upper truncated is 4100 and the lower truncated is 3900, and the standard deviation = 100.

#### *4.1. Impact of the p on SA\*, SAr\*, and F(T)*

In this subsection, the we investigate the impact of the returned products percentage, *p*, on the optimal stock capacity for new products, *SA*\*, reconditioned products, *SAr*\*, and the total profit, *F*(*T*). Thus, we vary the value of *p*, and by using the optimization algorithm, we determine the corresponding *SA*\*, *SAr*\*, and *F*(*T*). The simulations' results are illustrated in Table 2. In order to improve the visualization of the obtained results in the analyzed part, we add two Figures, Figures 4 and 5, to the presented Table.

**Table 2.** Optimal decision variables in function of *p*.


**Figure 4.** *SA*\* and *F*(*T*) in function of *p*.

**Figure 5.** *SAr*\* and *F*(*T*) in function of *p*.

As shown in Table 2, the total profit is, at the first time, proportional to the percentage of returned products. Then, as the percentage of returned products increases, the profit again decreases and leads to losses. Indeed, it is observed that we have an optimum when *p* is equal to 50%. This is explained by the fact that when *p* increases, the number of worn products returned, *R*(*t*), increases and replenishes the recovery inventory. Thus, when *p* is below 50%, the value of profit function increases, which is normal, and the quantity of returned products increases and replenishes the recovery inventory, then the reconditioned products are satisfied. Nevertheless, when *p* it is above 50%, many of the returned products will be stored in the recovery inventory, which leads to overstocking and generates a high storage cost, and thus explains the losses. The optimal manufacturing stock, *SA*\*, remains constant, which means that it does not depend on *p* (see Figure 4). In fact, the production of new parts is independent from the production of reconditioned parts. The optimal reconditioning stock, *SAr*\*, increases when *p* increases (see Figure 5), and this is explained by the need of the worn parts for supplying the demand of reconditioned parts in Market 2. If more worn parts are returned to the production system, more reconditioned parts will be available in the stock, *SAr*, to satisfy the demand, *dAr*(*t*). This study allows a firm leader to manage, in an optimal way, the manufacturing, reconditioning, and stocking when the quantity of returned worn products varies. Furthermore, when the leader has the option to set the quantity of the returned worn products, he can improve the system management by proposing the optimal quantity of the returned worn products that maximize the profit. Therefore, in the next studies, we add *p* to the decision variables.

#### *4.2. Impact of the W on* SA*\*, SAr\*, p\*, and F(T)*

In this subsection, we investigate the impact of the length of the mandatory carbon period, *W*, on the optimal stock capacity for new products, *SA*\*, reconditioned products, *SAr*\*, returned products percentage, *p*, and total profit, *F*(*T*). Thus, we vary the value of *W*, and by using the optimization algorithm, we determine the corresponding *SA*\*, *SAr*\*, *p*\*, and *F*(*T*). The simulations' results are illustrated in Table 3. In order to improve the visualization of the obtained results in the analyzed part, we add the Figure 6 to the presented Table.


**Table 3.** Impact of length of the mandatory carbon period *W* on *SA*\*, *SAr*\*, *p*\*, and *F*(*T*).

**Figure 6.** *SAr*\* and *SA*\* in function of *W*.

As it is shown in Table 3, the larger the mandatory carbon period is, the more the total profit decreases, which is explained by lower production in a large period. Also, the optimal of the length of the period of carbon is given (the bold line). Indeed, when the mandatory carbon period is low, the limitation for producing either for new or reconditioned products is low, thus the production is less limited, and the number of satisfying products increases, and then the lost sales costs decreases (see Figure 6). Consequently, to hold the produced products, the optimal storage capacities for both types of products increase. Thus, the demands will be satisfied more, and of course, the profit increases. This study allows a company manager to propose optimal decisions on the manufacturing, reconditioning, stocking, and quantity of returned worn products when a mandatory carbon period is imposed. Indeed, the government imposes a mandatory carbon period when the pollution reaches a level that is considered harmful to humans. In this period, the manager has to respect a fixed limit of carbon emissions. Consequently, the production is limited, and the manager has to find the optimal planning of production, stocking, and returning of worn products.

#### *4.3. Impact of the q<sup>l</sup> on SA\*, SAr\*, p\*, and F(T)*

In this subsection, the we investigate the impact of the limited quantity of carbon emission, *q<sup>l</sup>* , on the optimal stock capacity for new products, *SA*\*, reconditioned products, *SAr*\*, returned products percentage, *p*, and the total profit, *F*(*T*). Thus, we vary the value of *q<sup>l</sup>* , and by using the optimization

algorithm, we determine the corresponding *SA*\*, *SAr*\*, *p*\*, and *F*(*T*). The simulations' results are illustrated in Table 4. In order to improve the visualization of the obtained results in the analyzed part, we add the Figure 7 to the presented Table.


**Table 4.** Impact of the limited quantity of carbon emission, *q<sup>l</sup>* , on *SA*\*, *SAr*\*, *p*\*, and *F*(*T*).

**Figure 7.** *SAr*\* and *SA*\* in function of *q<sup>l</sup>* .

As it observed, when the limited quantity of carbon increases, the profit slightly increases, but both *SA*\* and *SAr*\* decrease (see Figure 7). In fact, when the limited quantity of carbon emission is high, the production either for the new or reconditioned products are less limited, and then both types of products are more available, thus, both demands *dA*(*t*) and *dAr*(*t*) are better satisfied. Moreover, when both types of products are available, the stocking decreases as the demands are satisfied. However, when the limited quantity of carbon emission is low, the production either for the new or reconditioned products are very limited. Thus, both types of products are less available and then it shall increase the storing to satisfy the demands. This study allows a firm manager to propose optimal decisions on the manufacturing, reconditioning, stocking, and quantity of returned worn products when the government imposes a limited quantity of carbon emission.

#### *4.4. Impact of the ct on SA\*, SAr\*, p\*, and F(T)*

In this subsection, we investigate the impact of the cost of carbon emission, *ct*, on the optimal stock capacity for new products, *SA*\*, reconditioned products, *SAr*\*, returned products percentage, *p*, and the total profit, *F*(*T*). Thus, we vary the value of *ct*, and by using the optimization algorithm, we determine the corresponding *SA*\*, *SAr*\*, *p*\*, and *F*(*T*). To study the influence of the cost of carbon emission, we changed the parameters' values for *cs<sup>A</sup>* = 0.00002, *csAr* = 0.00002, and *cr* = 0.0000001 monetary units. The simulations' results are illustrated in Table 5. In order to improve the visualization of the obtained results in the analyzed part, we add the Figure 8 to the presented Table.

**Table 5.** Impact of the cost of carbon emission on *SA*\*, *SAr*\*, *p*\*, and *F*(*T*).


**Figure 8.** *SAr*\* and *SA*\* in function of *ct*.

The first observation in Table 5 is that the total profit obviously decreases when the cost of carbon increases. Indeed, when the cost of carbon emission increases, the system produces less new parts in order to decrease the cost of emissions, and then, optimal stock, *SA*\*, decreases (see Figure 8). On the other hand, in order to decrease the cost of emissions, the system produces more reconditioned parts instead of new ones. That explain the increase of optimal stock, *SAr*\*, when the cost of carbon emission increases (see Figure 8). This study allows a firm manager to propose optimal decisions on the manufacturing, reconditioning, stocking, and quantity of returned worn products when the government increases the carbon emission cost.

#### **5. Conclusions**

This paper has proposed a new design of a manufacturing and reconditioning system, taking into account two carbon regulations: carbon tax and periodic mandatory emission. New and reconditioned products are distinguished and sold in different markets. Our work provides, in the framework of this system, a new and quick strategy to quantify the production of parts under the periods of carbon limitation. To make this green manufacturing more realistic, we have considered the failure and repair of the machines, and the amount of returned products (worn products) depends on the previous sales. The demands are stochastic. To simulate the proposed system and to faithfully describe the system behavior, we have developed a mathematical model based on a discrete flow model. Then, an optimization method based on an evolutionary algorithm was developed, to determine the optimal quantity of the returned worn products, and the stock capacity for new and reconditioned products. Four numerical results were obtained by using the optimization algorithm, which help a firm manager to make decisions on the production and storing of new and reconditioned parts. The results have provided optimal strategies of production, storing, and recovery of worn parts in function of the mandatory carbon period length, emission limit, and carbon cost.

In future research, we will consider the mentioned limitations and the carbon emissions exceeded by transportation from stores to markets.

#### **6. Discussion**

This work concerned the consideration of a periodic mandatory carbon policy that is combined with carbon tax. A new approach was developed to determine the quantity of new and reconditioned products to produce when the government imposes a limited carbon emission in a random period, the optimal levels of the manufactured-reconditioned stock, and the optimal percentage of returned worn products that maximizes the total profit. Furthermore, the proposed approach combined different options that make the study more realistic, such as: distinction between new and reconditioned parts, the dependence of worn products on the sales in the previous periods, and the variability of the machines repair periods. An efficient optimization method based on evolutionary algorithms was developed that returns the optimal values of the vector, which maximizes the total profit function. The studies provided in this work allow a company manager to propose optimal decisions on the manufacturing, reconditioning, stocking, and quantity of returned worn products when a mandatory carbon period is imposed or when the government increases the carbon emission cost. This work is not exempt from limitations. First, it is considered that the reconditioning process is assured by one machine, but in practice, usually it consists of several machines that are subject to random failures. Second, only one type of product was considered, whereas in practice, almost all companies are producing multiple products. Therefore, in our future research, we will extend the model by assuming that the manufacturing process as well as the reconditioning process consists of several activities carried out at several workplaces. Also, we will consider a multi-product system of production and reconditioning.

**Author Contributions:** S.T. is the responsible for conducting this research. S.T. proposed the main paper contribution. S.S. developed the mathematical model. C.S. developed the optimization method and algorithms. N.S. supervised this work and contributed to revise this paper. 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.

#### **Appendix A**

**Algorithm 1**: Function near.

```
Near (M[n], x, S[n], y)
01: M[n], x, S[n], y
02: Int i, result; Double D;
03: D = |TT [0] − x| + |SS [0] − y|
04: Result = 0;
05: For (i = 1; i < n; i++){
06: If (|TT[i] − x| + |SS[i] − y|; <= D)
07: D= |TT[i] − x|+|SS[i] − y|
08: Result = i;
09: End if
10: End for
11: Return Result
```
#### **Appendix B**

**Algorithm 2**: Production planning under mandatory carbon emission.

**Production plan** (*uA*(*t*), *uAr*(*t*)) 01: *k* = *0*; 02: **If** (*uA*(*t*)! 0 and *uAr*(*t*)! = 0) 03: PSC = *uA*(*t*)/(*uA*(*t*) + *uAr*(*t*)) 04: PAC = *uAr*(*t*)/(*uA*(*t*) + *uAr*(*t*)) 05: **For** (*i* = 1 to *uA*(*t*)) 06: **For** (*j* = 1 to *uAr*(*t*)) 07: **If** (*i*·*qpm* + *j*·*qpr* < = *qL*) 08: *Pm*2*co* = *i*/(*i* + *j*) 09 *PR*2*CO* = *j*/(*i* + *j*) 10: *h*[*k*] = *Pm*2*co* 11: *B*[*k*] = *PR*2*CO* 12: *v*[*k*] = *i* 13: *c*[*k*] = *j* 14: *k*++ 15: **End if** 16: **End for** 17: **End for** 18: *lignenear* = **Near** (*h*, *PSC*,*B*, *PAC*) 19: *i* = *v*[*lignenear*], *j* = *c*[*lignenear*] 20: **End if** 21: *uA*(*t*) = *i* 22: *uAr*(*t*) = *j*

#### **Appendix C**

**Algorithm 3**: Pseudo-code of couple of coordinates, iterative neighborhood.



## **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Spray and Aerosolised pH-Neutral Electrochemically Activated Solution Reduces** *Salmonella* **Enteritidis and Total Bacterial Load on Egg Surface**

**Sangay Tenzin <sup>1</sup> , Sergio Ferro <sup>2</sup> , Samiullah Khan <sup>3</sup> , Permal Deo 4,\* and Darren J. Trott 1,\***

	- samiullah.khan@adelaide.edu.au

**Abstract:** The effectiveness of sprayed and aerosolised pH-neutral electrochemically activated solutions (ECAS) containing 150 mg/L of free available chlorine in reducing total bacteria load and artificially inoculated *Salmonella enterica* serotype Enteritidis 11RX on eggs surfaces was investigated. Treatment groups included untreated control, sodium hypochlorite (positive control), sprayed and aerosolised water and sprayed and aerosolised ECAS. Sprayed ECAS (150 mg/L, 45 s) showed a significant reduction in total bacterial load (2.2 log reduction, *p* < 0.0001) and *S.* Enteritidis (5.4 log reduction, *p* < 0.0001) when compared with the untreated control. Aerosolised ECAS (120 s) was effective in reducing both the total bacterial load (1.4 log reduction, *p* < 0.01) and *S.* Enteritidis (4.2 log reduction, *p* = 0.0022). However, aerosolised ECAS (60 s) only significantly reduced *S.* Enteritidis counts (2.8 log reduction, *p* < 0.0008), indicating that a longer time for bacterial reduction during fogging sanitisation is needed. Tests performed with one egg per oscillating tray were more effective in reducing both the total bacterial load and the *S.* Enteritidis counts than those with three eggs per oscillating tray. Sprayed ECAS (45 s) and aerosolised ECAS (120 s) did not deteriorate the egg cuticle integrity (∆*E*ab \* ), which was evaluated using Cuticle Blue dye solution and colour intensity measurement. Overall, both the reduction in total bacteria counts and *S.* Enteritidis from the egg surface and retention of cuticle integrity suggest that sprayed and aerosolised ECAS could be used as alternative sanitising approaches to improve the food safety aspect of table eggs.

**Keywords:** pH-neutral electrochemically activated solution; total bacterial count; *Salmonella* Enteritidis; egg cuticle integrity

#### **1. Introduction**

Pathogenic serotypes of *Salmonella* are a major cause of foodborne diseases worldwide. The annual proportion of food origin salmonellosis in Australia is about 40,000 out of an estimated total of about 4.1 million foodborne gastroenteritis cases [1]. *Salmonella-* associated foodborne illnesses have risen during the past 20 years and the rate of salmonellosis in Australia is much higher compared to economically similar countries [2]. It has been estimated that foodborne illnesses due to *Salmonella* spp. have caused up to 35% of hospitalisations and 28% of mortalities [1], and the hospitalisation and death cases were higher in comparison to other foodborne illnesses [1]. Among the salmonellosis cases of foodborne origin, raw eggs and egg-based products have the highest frequency [3–5]. For example, between 2001 and 2016, 50% of all foodborne *Salmonella* illnesses in Australia were attributed to the consumption of contaminated eggs [6]; 84% of these cases were caused by *Salmonella enterica* subspecies *enterica* serotype Typhimurium (*S.* Typhimurium) and three cases were

**Citation:** Tenzin, S.; Ferro, S.; Khan, S.; Deo, P.; Trott, D.J. Spray and Aerosolised pH-Neutral Electrochemically Activated Solution Reduces *Salmonella* Enteritidis and Total Bacterial Load on Egg Surface. *Appl. Sci.* **2021**, *11*, 732. https://doi. org/10.3390/app11020732

Received: 21 December 2020 Accepted: 11 January 2021 Published: 13 January 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/).

caused by *Salmonella enterica* subspecies *enterica* serotype Enteritidis (*S.* Enteritidis) [6]. In other countries, salmonellosis is caused predominantly by the serotype Enteritidis [7–9]. As egg consumption in Australia is approx. 245 eggs per capita and growing [10], the industry is continuously exploring alternative means to address *Salmonella* contamination.

Unlike the one-off input costs for the establishment of farm infrastructure and human resources, farms and industries incur ongoing costs for egg hygiene and egg safety management. Washing and disinfection of egg surface are the key steps involved in egg production to reduce the risk of egg-related foodborne illnesses and to maintain consumer confidence on the microbiological safety of eggs. Elimination of pathogenic bacteria from the egg surface is achieved using several techniques and many eggshell sanitisation methods are also employed to reduce the contamination of eggs by *Salmonella* in commercial egg production premises.

Protocols employed for the reduction of *Salmonella* can be broadly classified into thermal and non-thermal disinfection procedures. Thermal disinfection, such as egg pasteurisation, is a highly effective method, but negatively affects egg proteins and rheological properties [11]. The most common non-thermal processes employ quaternary ammonium compounds (QAC) and chlorine-based chemicals [12] to sanitise eggs after washing with a high pH (11.0) detergent solution at a temperature above 40 ◦C. Unfortunately, bacteria may develop resistance to QAC [13], which in turn induces selection of genes for co-resistance to other antimicrobials [14], thus limiting its use. In the case of chlorine-based sanitisers, besides the development of bacterial resistance due to its persistent usage [15], the accumulation of organic load from dirt, manure and broken eggs reduces the chlorine concentration, affecting the efficacy of chlorine-based sanitisation. Moreover, due to the environmental impact caused by chlorine-based by-products and the problems with wastewater disposal, its usage in the food industry is limited. Other decontamination methods used are ultraviolet (UV) irradiation of eggs after washing and formaldehyde fumigation. However, the antibacterial activity of UV irradiation protocol is limited to the egg surface directly exposed to UV rays [16], while formaldehyde is a known occupational health hazard and a carcinogenic product [17].

Since occupational health safety and environmental regulations continually push towards safer eggshell sanitisers, electrochemically activated solutions (ECAS) (also called electrolysed oxidizing (EO) water) could be a potential alternative for eggshell cleaning and disinfection. The three forms of ECAS (acidic, slightly acidic and neutral) have been previously assessed for the sanitisation of table eggs in safe food production [18–21] and fertile eggs for quality production of chicks [22]. In most of the available research, a twostep process for ECAS disinfection of eggs was followed. In the initial washing step, dirt and debris are washed off the egg surface with water or alkaline detergent, followed by the ECAS disinfection. Bialka et al. [19] reported that immersion washing of eggs with acidic ECAS significantly reduced *S.* Enteritidis and *Escherichia coli* from the egg surface but also damaged the egg cuticle layer. On the other hand, the spray washing of eggs with slightly acidic EO water reduced total aerobic bacteria without negatively affecting the cuticle layer [22]. In other studies, an immersion washing with pH-neutral ECAS was not effective in reducing the total bacterial load [23], while a spray wash significantly decreased the level of *Listeria monocytogenes* without affecting the egg cuticle layer [20]. More recently, Medina-Gudino et al. [24] reported that pH-neutral EO spray treatment for 30 s significantly reduced *S.* Enteritidis and *E. coli* loads on the egg surface without adversely affecting egg cuticle integrity.

In this study, we explored the potential of pH-neutral ECAS (150 mg/L of free available chlorine (FAC)) as spray and aerosol fog for the sanitisation of unwashed, visibly clean eggs, assessing the reduction in total bacterial counts and *Salmonella* Enteritidis, and its effects on the cuticle layer. ECAS with neutral pH still contains hypochlorous acid (HOCl) as the main oxidising component (active chlorine compounds also include hypochlorite ions and dissolved gaseous chlorine) [25–27] but is less corrosive and more durable than the acidic and slightly acidic forms.

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

#### *2.1. Chicken Egg Source and Selection*

Freshly laid eggs were sourced from Hy-Line Brown hens (aged between 36 and 40 weeks) raised in conventional battery cages housing individual hens at the School of Animal and Veterinary Sciences, The University of Adelaide. Eggs stored for 24 h at room temperature were visually screened for thermal cracks and dirt; relatively clean, intact and uniformly sized eggs were selected and randomly divided into groups to determine the effectiveness of sanitisation treatments on total microbial load and artificially inoculated *S*. Enteritidis load.

#### *2.2. ECAS Spray and Fog Disinfection Generation*

The electrochemically activated solution was sourced from Ecas4 Australia Pty Ltd. (Adelaide, Australia) and its physicochemical properties such as temperature, pH, and oxidation/reduction potential (ORP) were measured using a handheld meter (Model MC-80, TPS Pty Ltd., Brisbane, Australia). Free and total available chlorine was measured using a Free Chlorine Checker® HC-HI701 and a Total Chlorine Checker® HC-HI711 (Hanna Instruments, Melbourne, Australia), respectively. ECAS was stored at 4 ± 1 ◦C and used within one week of preparation.

The working concentration of ECAS (pH ≈ 6.8–7.0, 150 mg/L of free available chlorine (FAC)) was freshly prepared each time prior to sanitisation experiment. We have previously optimised the aerosolised ECAS at 150 mg/L of FAC to significantly reduce the total microbial load in an animal farm environment, and for this reason we have chosen this concentration [28]. For the spray wash, ECAS was sprayed for 45 s using a handheld bottle sprayer. For the fog sanitisation, ECAS fog was generated using an ultrasonic humidifier that generates droplets sized between 1 and 3 µm in diameter (Ultrasonic Humidifier HU-85, Contronics Engineering B.V., Sint-Oedenrode, The Netherlands) as previously optimised [28].

#### *2.3. Effectiveness of Spraying and Fogging on Total Bacterial Load Reduction on Eggshell Surface*

For total bacterial counts, intact and visibly clean eggs were selected for each of the treatment groups (three eggs per treatment): unwashed control, NaOCl spray (~200 mg/L of FAC, positive control), water spray (45 s), ECAS spray (150 mg/L of FAC; 45 s), water fog (60 and 120 s) and ECAS fog (150 mg/L of FAC; 60 and 120 s). The disinfection procedures were performed in a biosafety cabinet (BSC) and eggs were placed in separate compartments on an oscillating tray for uniform exposure during spray-wash and fogging treatment. Each treatment was independently repeated for at least two times. In addition, ECAS (150 mg/L of FAC) fog treatment was compared for one and three eggs per oscillating tray, for 60 s and 120 s, respectively.

After treatment, individual eggs were immediately placed into a sterile Whirl-Pak bags (Nasco, Fort Atkinson, WI, USA) containing 5 mL of buffered peptone water (BPW), massaged gently (without breaking the egg) for 1 min, then the broth was transferred into 10 mL sterile tubes (SARSTEDT Australia Pty Ltd., Adelaide, Australia). Samples were centrifuged at 5444× *g* (MPW-351e Centrifuge, Med Instruments, Adelab Scientific, Adelaide, Australia) for 10 min, supernatant discarded, and the pellet was resuspended in 200 <sup>µ</sup>L of 1 <sup>×</sup> PBS. Aliquots (100 <sup>µ</sup>L) of 10-fold dilutions (up to 10−<sup>5</sup> ) of the samples were spread plated on plate count agar (PCA; CM0325, Oxoid, Melbourne, Australia) in duplicates and incubated overnight at 37 ◦C for the enumeration of colonies. Plates with 25 to 300 colonies were used for colony forming unit (CFU) calculation and data presented as log<sup>10</sup> CFU/egg.

#### *2.4. Salmonella Enteritidis Seeding on Outer Eggshell Surface*

#### 2.4.1. Pre-Wash of Eggs for *S*. Enteritidis Seeding

To understand the efficacy of ECAS on *Salmonella* load reduction, eggs selected for *Salmonella* seeding were washed as per wash steps specified in Gole et al. [29] before being inoculated with *S.* Enteritidis. Briefly, eggs were placed on an oscillating tray, which helped in exposing the entire eggshell surface, in a BSC and initially spray-washed with a 0.45% (*v/v*) solution of NaOCl (ThermoFisher, Melbourne, Australia; pH ≈ 12) at 40 ◦C for 45 s. Then, spray-sanitised with a 0.16% (*v/v*) solution of NaOCl at 32 ◦C for 22 s and left on the sterilised BSC bench to dry for 60 min. The eggs were sanitised to ensure the complete removal of the microbiota of the egg surface and to achieve a uniform *S.* Enteritidis colonisation of the egg surface.

#### 2.4.2. *S*. Enteritidis Preparation for Inoculation of Eggshells

*Salmonella* Enteritidis 11RX was used for this experiment. *S.* Enteritidis stored at −80 ◦C in 80% glycerol was plated on xylose lysine deoxycholate (XLD) agar (Oxoid CM0469) and incubated overnight at 37 ◦C to obtain isolated colonies. An inoculum was prepared by suspending colonies in phosphate buffered saline (1 × PBS) to obtain an absorbance (OD600 nm) value of 0.45. Viable *Salmonella* was enumerated by plating 10-fold serial dilutions of the inoculum on XLD agar and incubating overnight at 37 ◦C. After enumeration, a 200 mL inoculum containing ~10<sup>5</sup> CFU per mL was prepared in 1 <sup>×</sup> PBS.

For eggshell seeding, eggs were immersed for 90 s either in 1 × PBS (control) or in <sup>1</sup> <sup>×</sup> PBS containing ~10<sup>5</sup> CFU/mL of *S.* Enteritidis. Eggs were then placed into sterile zip lock bags and incubated at 37 ◦C. After 18–24 h post-inoculation, three eggs from each treatment were placed in separate Whirl-Pak bags containing 5 mL BPW and massaged for 1 min. Aliquots (100 µL) from a 10-fold serial dilution were spread plated on XLD and PCA media (in duplicate) and incubated overnight at 37 ◦C for enumeration of counts.

#### 2.4.3. Effectiveness of Spraying and Fogging on *S.* Enteritidis on Eggshell Surface

The *S.* Enteritidis-inoculated eggs were subjected to the same treatment as above (Section 2.3). A 10-fold serial dilutions were prepared for each treatment as above in 1 × PBS and aliquots were spread plated both on PCA and XLD agar media (in duplicates).

#### *2.5. Eggshell Cuticle Assessment*

Twelve eggs per each treatment (water and ECAS spray washing for 45 s and water and ECAS fogging for 2 min) were screened and selected based on colour intensity measured using a MiniScan EZ colourimeter (4500 L Spectrophotometer, Hunter Associates Laboratory, Inc., Reston, VA, USA). The selected eggs were treated as in Section 2.3 and dried in biosafety cabinet for 60 min. The eggs were stained with MST Cuticle Blue dye (MS Technologies, Kettering, UK) as described by Khan et al. [30] and the cuticle coverage was assessed using the ∆*E*ab \* method. The average of four readings of the *L\** (lightness), *a\** (red/green value) and *b\** (yellow/blue value) values, before and after staining, were used for the calculation of ∆*E*ab \* with Equation (1). A higher value of ∆*E*ab\* denotes a higher cuticle staining:

$$
\Delta E\_{\rm ab} \stackrel{\*}{=} \sqrt{[(\Delta L^\*)^2 + (\Delta a^\*)^2 + (\Delta b^\*)^2]} \tag{1}
$$

#### *2.6. Statistical Analysis*

Total bacterial and *S*. Enteritidis counts were expressed as log10/egg and the data were analysed using ANOVA in GraphPad Prism v.8 (GraphPad Software, San Diego, CA, USA). Since the bacterial counts were log transformed (subjected to normal distribution) and each experiment had an equal number of samples, a one-way ANOVA was performed to compare differences of means among untreated control versus different sanitising treatments followed by a Tukey's multiple comparison test. A *p* value of < 0.05 was considered statistically significant.

#### **3. Results**

#### *3.1. ECAS Spray and Aerosol Treatment Reduced Total Bacterial Count on Egg Surface*

The mean bacterial load (log<sup>10</sup> CFU/egg) for untreated and after sanitising treatments are presented in Table 1 and Figure 1. Eggs treated with sprayed water (45 s) (2.4 ± 0.1 log<sup>10</sup>

CFU/egg) showed no significant reduction (*p* = 0.3662) in total bacterial load when compared with the untreated control (2.2 ± 0.2 log<sup>10</sup> CFU/egg; Table 1). Significant reduction in total bacterial load was observed when the eggs were treated with sprayed ECAS (45 s) (2.2 log reduction, *p* < 0.0001) and aerosolised ECAS (120 s) (1.4 log reduction, *p* = 0.01) compared to untreated control (Table 1). Aerosolised water (60 s and 120 s) and aerosolised ECAS (60 s) showed no significant reduction in total bacterial load when compared with sprayed water (45 s) or the untreated control (Table 1). In addition, the spray versus aerosol techniques were compared for their effectiveness in reducing the total bacterial load (Figure 1A). Sprayed ECAS (45 s) showed a significant reduction (*p* < 0.0001) in total bacterial load when compared with sprayed water (45 s) (Figure 1A). Aerosolised ECAS (120 s) also showed a significant reduction (*p* < 0.01) in total bacterial load when compared with aerosolised water (120 s); however, no significant reduction was shown for the 60 s treatment group (Figure 1A).



Total bacterial counts were calculated as log<sup>10</sup> CFU/egg and presented as mean <sup>±</sup> standard deviation (SD); & log<sup>10</sup> reduction = (log<sup>10</sup> counts of untreated control)—(log<sup>10</sup> counts after sanitising treatment); a '+' sign means log increase in counts.

✱✱

✱✱

✱✱✱✱

**Figure 1.** *Cont*.

✱✱✱✱

**2021**, , x FOR PEER REVIEW 6 of 12

✱✱

✱✱✱✱

**Figure 1.** Effectiveness of sprayed and aerosolised ECAS on reduction of total bacterial load on egg surface. (**A**) Total bacterial load after the various sanitisation treatments. (**B**) Aerosolised ECAS sanitisation of individual (one egg/tray) and simultaneous (three eggs/tray) treated eggs for 60 s and 120 s. NaOCl—sodium hypochlorite (~200 mg/L FAC); ECAS—electrochemically activated solution (150 mg/L of FAC); ns—not significant; \*\* *p* < 0.01, \*\*\*\* *p* < 0.0001.

The effectiveness of aerosolised ECAS was further assessed when the eggs were sanitised simultaneously (three eggs/tray) or individually (one egg/tray) for 60 s and 120 s, respectively. Aerosolised ECAS (120 s) significantly reduced the total bacteria load for both one egg/tray (*p* < 0.0001) and three eggs/tray (*p* < 0.01; Figure 1B). The treatment of three eggs/tray with aerosolised ECAS (120 s) did not eliminate the total bacterial load completely; however, the level was not significantly different from the treatment of one egg/tray aerosolised ECAS (120 s) (Figure 1B).

#### *3.2. ECAS Spray and Fog Treatment Reduced S. Enteritidis on the Egg Surface*

*S.* Enteritidis counts (log<sup>10</sup> CFU/egg) are presented in Table 2 and Figure 2. Sprayed water (45 s) showed a significant reduction (1.0 log reduction, *p* = 0.0005; Table 2) when compared with the untreated control. All the ECAS treatments significantly reduced the inoculated *S.* Enteritidis counts on the egg surfaces compared to the untreated control: sprayed ECAS (45 s) (5.4 log reduction, *p* < 0.0001; Table 2), aerosolised ECAS (60 s) (2.8 log reduction, *p* = 0.0008; Table 2) and aerosolised ECAS (120 s) (4.2 log reduction, *p* = 0.0022; Table 2). A significant reduction in *S.* Enteritidis counts (1.1 log reduction, *p* < 0.0001) was also achieved with aerosolised water (120 s), whereas no significant reduction was observed with aerosolised water (60 s) when compared with the untreated control (Table 2). For the spray versus aerosol techniques comparison, sprayed ECAS (45 s) showed a significant reduction (*p* < 0.0001) in *S.* Enteritidis load when compared with the sprayed water (45 s) (Figure 2A). Significant reduction in *S.* Enteritidis load was observed for aerosolised ECAS (60 s, *p* < 0.01) and aerosolised ECAS (120 s, *p* < 0.01) when compared with respective aerosolised water treatments (Figure 2A). *S.* Enteritidis level was significantly reduced in spray water (45 s, *p* = 0.0005) and aerosolised water treatment groups (120 s, *p* < 0.0001), respectively.

**2021**, , x FOR PEER REVIEW 8 of 12

**Figure 2.** Effectiveness of sprayed and aerosolised ECAS on reduction of *S.* Enteritidis counts. (**A**) *S.* Enteritidis counts after the various sanitisation treatment. (**B**) Aerosolised ECAS sanitisation of individual (one egg/tray) and simultaneous (three eggs/tray) treated eggs for 60 s and 120 s. NaOCl—sodium hypochlorite (~200 mg/L FAC); ECAS—electrochemically activated solution (150 mg/L of FAC); ns—not significant; \*\* *p* < 0.01, \*\*\* *p* < 0.001, \*\*\*\* *p* < 0.0001.

Δ


**Table 2.** Effect of sanitisation treatments on the reduction level of *S.* Enteritidis counts on egg surface.

*S.* Enteritidis counts calculated as log<sup>10</sup> CFU/egg and presented as mean <sup>±</sup> standard deviation (SD); & log<sup>10</sup> reduction = (log<sup>10</sup> counts of untreated control)—(log<sup>10</sup> counts after sanitising treatment).

Aerosolised ECAS (120 s) significantly reduced *S.* Enteritidis counts for one egg/tray (*p* < 0.0001); however, no significant reduction trend was observed in the 3 eggs/tray treatment group (Figure 1B).

#### *3.3. ECAS Spray and Aerosol Treatments Did Not Affect Egg Cuticle Integrity* **2021**, , x FOR PEER REVIEW 9 of 12

The effect of sprayed and aerosolised ECAS treatments on cuticle coverage (∆*E*ab \* ) was measured (Figure 3). Sprayed ECAS (45 s) and aerosolised ECAS (120 s) did not show any significant difference when compared with the respective water controls.

Δ **Figure 3.** Effect of sprayed and aerosolised ECAS treatments on cuticle integrity (∆*E*ab \* ). Data presented as mean ± SD, n = 12; ns—not significant.

#### **4. Discussion**

This study assessed the effectiveness of sprayed and aerosolised ECAS in reducing total bacterial load and inoculated *S.* Enteritidis on egg surface and their effects on the cuticle layer. The bacterial load on eggshells is usually acquired through contamination from the farm environment; therefore, farm type and poultry housing system influence the egg surface total bacterial count [31–34]. Eggshells from conventional-caged hens usually harbour lower total bacterial counts compared to other housing systems [35]. The total bacterial load of 2.2 ± 0.2 log<sup>10</sup> CFU/egg found in this study was lower than those observed by Wall et al. [36] (2.7 log10) and Alvarez-Fernandez et al. [31] (2.3 log10) on eggshells from conventional-caged hens, probably because the eggs used in this study were laid by hens housed individually with a low density of hens in the shed (49 hens/shed).

ECAS spray-wash of eggs for 45 s and aerosolisation of individual eggs with ECAS for 120 s completely reduced the native total bacterial load and *S.* Enteritidis on eggshell surface. These two sanitisation approaches showed an effectiveness like that of a sodium hypochlorite (~200 mg/L of FAC) spray-wash for 45 s, which was used as a positive control. Although relevant literature is not available on the sanitisation of eggshell surface with ECAS fog, previous research has shown a significant reduction in total bacterial counts from the eggshell surface when the latter is spray-washed with pH-neutral EO water or acidic EO water [22,37]. In the present study, the effectiveness of aerosolised ECAS with one and three eggs per oscillating tray for 60 s and 120 s was also tested. ECAS fogging of multiple eggs (three simultaneously) as well as of a single egg for 120 s showed a significant reduction in total bacterial load when compared with the respective 60 s treatments, confirming that a longer treatment time may be more appropriate for total bacteria load reductions. It is noteworthy that a higher variation was observed in treating three eggs simultaneously, which could be attributed to an uneven distribution of the fog, as fogging was performed in a biological safety cabinet that quickly sucked out the surrounding air including the fog.

Previous studies have assessed the effectiveness of various forms of electrolysed water (acidic, slightly acidic, neutral, alkaline) by immersion or spray-washing of eggshell surface inoculated with *S.* Enteritidis [24,38–40]. In this study, we assessed a pH-neutral EO water in the form of fog in addition to spray-washing. The significant microbial reduction observed in our study for sprayed water and aerosolised water (120 s) could be due to the dislodgment of the inoculated *S.* Enteritidis from the egg surface; however, these techniques were ineffective in the reduction of total bacterial load. Sprayed ECAS (45 s) and aerosolised ECAS (60 s and 120 s) showed significant reduction in *S.* Enteritidis counts. Medina-Gudino et al. [24] reported that pH-neutral EO spray treatment (30 s) significantly reduced *S.* Enteritidis (>1.45 log<sup>10</sup> CFU/egg) on egg surface. In another study [39], a slightly acidic (pH 6.53) EO water (15 mg/L of FAC) used at 20 and 45 ◦C showed a 4.2 log<sup>10</sup> CFU/mL reduction of *S.* Enteritidis. Aerosolised ECAS (120 s) showed a significant log reduction of *S.* Enteritidis for one egg/tray treatment but no significant reduction was observed in the three eggs/tray treatment group. As previously mentioned, the higher variation in the three eggs/tray treatment could be due to the fog being quickly sucked by the BSC and not moistening the egg surface with enough FAC to reduce the bacterial cells.

The main purpose of egg washing is to reduce the bacterial load on the egg surface; however, one of the major concerns in the egg washing process is the likely damage of the cuticle layer. In the present study, no difference was observed in ∆*E*ab \* values in eggs sprayed and aerosolised with ECAS, indicating that cuticle integrity was not altered. In a recent study, Medina-Gudino et al. [24] reported no differences in the overall ∆*E*ab \* values when eggs were treated with pH-neutral EO. In contrast, for eggs washed with alkaline and acidic EO, changes in the *a\** and *b\** values were observed, hence the reduction in overall ∆*E*ab \* values, indicating damage to the egg cuticle [19].

Although electrolysed water is an environmentally friendly, non-hazardous sanitiser with proven antibacterial efficacy against foodborne pathogens such as *E. coli* O157: H7, *L. monocytogenes* [41] and *S.* Enteritidis on shell eggs [38,39], it is not currently being considered in commercial decontamination settings to obtain pathogen-free eggs. The reasons for this are probably related to the corrosiveness of acidic formulations towards steel surfaces and the lack of knowledge of consumers on the impact of chemicals for disinfection of environments.

The current scenario shows a shift in table egg production towards free-range systems [42]. For instance, in Australia, free-range production has increased from 39% in 2015 to 45% in 2018 [10], driven by consumer demand for bird welfare and access to a range area. However, this approach poses risks to public health as the eggshell bacterial load, including counts of *Campylobacter* and *Salmonella* spp. [42], in free range production systems is considerably higher than in conventional battery cage systems [42,43]. The problem is further aggravated by the increase in the consumption of meals consisting of

raw egg products [44]. Therefore, egg producers should consider eggshell cleaning and disinfection of table eggs as a priority to produce safe eggs and maintain consumer confidence. Besides issues with wastewater disposal, commonly used chlorine-based egg washing requires intensive monitoring of water temperature, pH, and chlorine concentration to retain its optimal efficacy and eggshell cuticle integrity. This suggests that electrolysed water could fill an important niche if washing and disinfection systems can be developed on a commercial level.

#### **5. Conclusions**

pH-neutral ECAS (150 mg/L of FAC), when used in the form of a spray (45 s) or as an aerosol (120 s), allows significant reductions in total bacterial load and *S.* Enteritidis counts from contaminated egg surface while retaining cuticle integrity. However, aerosolised ECAS (60 s) did not show significant reduction in either the total bacterial load or *S.* Enteritidis counts, suggesting a longer fogging time is needed to sanitise the eggs. The aerosolised ECAS sanitisation technique could be incorporated into cleaning and disinfection protocol to improve egg safety without the use of hazardous biocidal agents. Moreover, this disinfection protocol is easily implementable as ECAS can be easily generated on site with automated controls for FAC concentration and pH measurements. For this process to lead to large-scale tests and industrial implementation, additional testing needs to be performed, including the elimination of other pathogenic bacteria from eggshells, interior egg quality and consumer sensory evaluation, as per the regulatory requirements.

**Author Contributions:** Conceptualization, S.T., P.D., D.J.T.; Methodology, S.T., S.K., P.D.; Formal analysis, S.T., P.D.; Resources, S.F., S.K., P.D., D.J.T.; Data curation, S.T.; Writing—original draft preparation, S.T., P.D.; Writing—review and editing, S.T., S.F., S.K., P.D., D.J.T.; Visualization, S.T.; Supervision, P.D., D.J.T.; Funding acquisition, D.J.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** S.T. was supported by the Endeavour Postgraduate Scholarship.

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

#### **References**


## *Article* **Automated Chlorine Dosage in a Simulated Drinking Water Treatment Plant: A Real Case Study**

## **Javier Gámiz 1,2 , Antoni Grau 1,\* , Herminio Martínez <sup>3</sup> and Yolanda Bolea <sup>1</sup>**


Received: 21 May 2020; Accepted: 8 June 2020; Published: 11 June 2020

## **Featured Application: This research was applied in the Sant Joan Despí drinking water treatment plant in Barcelona city, providing drinking water to 1,900,000 inhabitants. The tool was validated with real data obtained in the plant in di**ff**erent conditions of e**ffl**uents.**

**Abstract:** In this paper, we present a simulator of a drinking water treatment plant. The model of the plant was based in hydraulic and matter transportation models. In order to not introduce more inaccuracies in the simulation, the control system was based in the real equipment deployed in the plant. This fact was the challenging part of the simulator, and an accurate design is presented in this research, wherein the sampling time had to be limited to interchange data between the SCADA in the plant and the simulator in real time. Due to the impossibility to stop the plant when testing the new control strategy, a simulator implemented the plant behavior under different extreme conditions. The validation of the simulator was performed with real data obtained from the plant.

**Keywords:** drinking water treatment plant; advection–diffusion–reaction; chlorine dosage; simulator validation; industrial application

#### **1. Introduction**

Drinking water is a necessary but scarce good. Despite the fact that water is the most abundant and common substance on our planet—as it covers 70% of its surface—96.5% of it is contained in the oceans. Of the remaining 3.5%, approximately 2.5% is found in the polar ice caps and glaciers and only 0.61% is liquid fresh water. Of the latter, around 0.98% is found in underground aquifers, which are difficult to access, while only 0.009% constitutes fresh surface water (rivers and lakes). Furthermore, only 0.003% of the total is fresh water available to be used for residential purposes. That is, if the Earth's total water were a 100 L container, only half a teaspoon of water would be suitable for human consumption [1]. The sources of water pollution can be natural (rain, decomposing vegetable matter, soil erosion, etc.) or anthropogenic (activity livestock, by-products of industrial activity, home waters, etc.), but both give rise to water that does not meet the necessary requirements to ensure its potability. The basic water treatment processes include several stages: coagulation, flocculation, particle separation (sedimentation/flotation), filtration, and disinfection (chlorination/ozonation) [2]. In many of these stages, the addition of chemicals into the flow of water to be treated is performed, and it is at this very point that this research provides a contribution with respect to correct dosage and control. Of all the treatments mentioned above, this paper focuses on disinfection. This process

attempts to destroy or inactivate pathogenic organisms present in water, mainly bacteria, viruses, and protozoa [3].

Disinfection treatments can be physical (gamma radiation, X-rays, ultraviolet radiation, thermal sterilization, etc.) or chemical (heavy metals, acids or bases, halogens, ozone, permanganate, etc.), with the latter being the most common form of treatment. Among the chemical reagents, chlorine and its derived compounds are the most widely used disinfecting agents worldwide. Many organisms regulate the residual chlorine values and they depend on the end use of the water. Thus, for drinking water, it is recommended that the residual free chlorine be between 0.2 and 1 ppm (parts per million), while in the case of swimming pools and spas, it should be kept between 1.5–3.0 ppm. However, these values are general and each competent body has determined its own thresholds.

The use of chlorine as a disinfecting agent began in the early 20th century and went on to complete the filtration process, which was already widely used. The most common chlorine family products for water disinfection are chlorine gas, chloramines, sodium hypochlorite, and calcium hypochlorite. Chlorine (Cl2) is a yellowish green, denser than air, toxic gas. It is a very oxidizing product that reacts with many compounds [4]. In fact, the most widely used chemical in the world is chlorine.

In the presence of humidity, chlorine is extremely corrosive and therefore the conduits and the materials in contact with it must be made of special alloys. Chlorine vapor is irritating by inhalation and can cause serious injury if exposed to high concentrations. Chlorine management must therefore be carried out by specialized staff, and very effective control and alarm systems are necessary. For these reasons, the use of hypochlorites in solution or in solid form is preferable. Sodium hypochlorite (NaClO) in solution is a disinfectant that has been used since the 18th century and is popularly known as bleach. At an industrial level, it is obtained by reacting chlorine gas with a sodium hydroxide solution. After the reaction, greenish yellow aqueous solutions are obtained, which have a determined concentration of active chlorine per liter [5]. It is marketed in solutions with concentrations between 3% and 15% by weight. Sodium hypochlorite is a very powerful and unstable oxidant, and thus a solution of 100 g of active chlorine per liter, after being stored for 3 months, can contain 90 g or even less. Calcium hypochlorite (Ca(ClO)2) is a white solid with a content of between 20% and 70% active chlorine. It is highly corrosive and can ignite on contact with certain acidic materials. However, it has two advantages over sodium hypochlorite: its higher chlorine content and its greater stability. To be used, it is diluted with water to obtain a more manageable concentration solution, for example 2%. When Cl<sup>2</sup> dissolves in water, it rapidly hydrolyses to generate hypochlorous acid and hydrochloric acid.

$$\text{Cl}\_2 + \text{H}\_2\text{O} \qquad \leftrightarrow \qquad \text{HClO} + \text{HCl}$$

In the case of hypochlorites, the dissociation of both salts occurs according to the following equations:

$$\begin{array}{cccc} \text{NaClO} + \text{H}\_{2}\text{O} & \leftrightarrow & \text{NaOH} + \text{HClO} \\ \text{Ca(ClO)}\_{2} + 2\text{H}\_{2}\text{O} & \leftrightarrow & \text{Ca(OH)}\_{2} + 2\text{HClO} \end{array}$$

Therefore, hypochlorous acid, which is actually the disinfectant variety, ends up forming in either case chlorine, sodium hypochlorite and calcium hypochlorite.One of the disadvantages of using chlorine and its derivatives is that it reacts with a large amount of organic matter and gives rise to trihalomethanes (THM), many of which have been shown to be toxic or carcinogenic [6]. Another drawback is the formation of chlorophenols in waters containing phenols, which would lead to bad odors. Chlorine also reacts with ammonia dissolved in water to form chloramines. These products also have some disinfecting power, although they are approximately 25 times less effective than free chlorine. However, their residence time in the water is long and, for this reason, they have sometimes been used as a reserve for residual chlorine. They present two major drawbacks: they can give rise to odors and flavors and are potentially toxic in a chronic way. Figure 1 shows an interesting representation of the risks associated with disinfection with chlorine. We adapted this figure from [7], indicating the optimal point for chlorination; due to the dynamics of chemical and microbiological risks, the optimal point of chlorination is at the intersection, in order to avoid both risks. It is clear that

not dosing can lead to a very high risk of infection (microbiological risk), but in any case, overdosing is not a valid solution, as it does not guarantee the elimination of health risks, since it favors the increase of chemical risks. An insufficient dose can cause vomiting or diarrhea (microbiological risk) and an overdose has carcinogenic effects (chemical risk); thus, the objective is an optimal dose. At very high levels of chlorine, the microbial risk increases, as taste and odor may cause the use of unsafe supplies.

**Figure 1.** Relationship between chlorination level and associated risk, adapted from [7].

On the basis of the foregoing, it can be deduced that chlorine (and derivatives), in addition to reacting with microorganisms, also reacts with other matter dissolved in the medium: organic matter, metals (iron, manganese), and mainly plastic derivatives. For this reason, to achieve a certain level of residual chlorine, the necessary amount to be added is much higher than the residual obtained. Therefore, before deciding on the dose of chlorine to be used to disinfect, the demand for chlorine must be determined, that is, the amount of chlorine consumed within the tank until the residual appears at the end of the plant, prior to entering the water distribution step.

Regarding the legislation of drinking water, there are some regulations at different levels. The higher level in Europe is dictated by the European Parliament and of the Council in its Directive 2008/105/EC on environmental quality standards in the field of water policy [8]. This legislation guarantees that water intended for human consumption can be consumed safely and without danger to health. Specifically, in Spain (valid in the current drinking water treatment plant (DWTP) location), the regulation is articulated through royal decrees RD 140/2003, RD 314/2016, and RD 902/2018 [9], specifying that the residual chlorine value must be between 0.2 and 1.0 ppm at all points in the supply network. In short, this legislation will be replaced by the future European Directive on Drinking Waters, which is at the moment in the final revision phase, with the date of possible approval at the end of 2020 and with mandatory entry in 2022. This new standard requires a series of new requirements and is mainly aimed at avoiding unnecessary water loss and helping to reduce the carbon footprint of the EU member states. The new water directive emerges to achieve the 2030 sustainable development goals (Goal 6) and the Paris agreement on climate change. Information will be key to increasing consumer confidence in tap water, and thus information on the quality and supply of drinking water in each area should be provided on the Internet. The idea is that the more information there is, the more confidence there will be about this resource that comes out of our taps, and therefore there will be a lower purchase of plastic bottled water, reducing the waste from this material.

The main goal of this research was to develop a simulator of a drinking water treatment plant (DWTP). New technologies deployed in this kind of plant allow for the automatization of the chlorination process, and the use of this novel simulator will allow for the testing of any chlorination strategy without involving such a critical infrastructure. The plant cannot stop delivering drinking water to a large portion of the population and the tests cannot be done with the plant at full capacity. The simulator will have the ability to be connected in real-time with automatic chlorination system

(specific hardware in the plant), and all this hardware equipment does not need to be simulated, thus enhancing the reliability of the whole process.

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

#### *2.1. Simulator Design*

This is one of the most important contributions of the presented research. The design of the simulator starts with defining the requirement, and proposing a novel architecture that will be able to fulfill them. This proposal involves the software implementation of different mathematical models (hydraulic and transportation models) and a hardware interconnection between plant computers, that is, programmable logic controller (PLC) and the main computer executing the simulator.

The chlorination process in the plant has been subjected to different phases of improvement at the automation level in recent years. The problems inherent in the automatic control have been diverse: variability in the chemistry of chlorine, punctual existence in the water to deal with significant amounts of ammonium, large dead times due to the hydraulics of its facilities, or disturbances caused by sudden changes in flow. The use of a simulator that reflects in the most realistic way the behavior of the chlorination process greatly facilitates the design of the control system in the effluent and allows decisions to be made on the optimization of the actual installation. Simulation has been seen over time as a modelling tool that has a very broad development and does not require sophisticated mathematics or statistics to develop a model and its usage [10]. On the other hand, the realization of experiments or tests on the real plant was expensive because of the need to have additional electronics and infrastructure, and moreover it mainly is especially dangerous because this is the last process of the DWTP and there is no correction mechanism in case of error (overdosage or underdosage of chlorine in tap water). The DWTP cannot be closed for testing new chlorination strategies, and the large duration of tests cannot be afforded by all the citizens and small industries that need drinking water uninterruptedly. Therefore, we decided to use a simulator that would allow for the study and analysis of the most appropriate online controller to later transfer the simulation results to the real control system of the plant.

When selecting a simulator to implement the desired model, we defined four strong requirements due to the complexity of the plant and the implications of stopping it:


In recent years, there is a growing interest by large automation companies to integrate computational fluid dynamics (CFD) software into online or real-time control scenarios. Through analyzing different software or platforms existing in the market to be used as a simulator in the chlorine automatic dosage, we finally decided to implement our own simulator, mainly for the simplicity of equations that we developed in this research. With our proposal, it is possible to reduce the sampling time to near 1 second. This fact is very important in order to have a good description and knowledge of water behavior and chlorine diffusion into tanks. After accurate research among many commercial simulators (18 software packages in total), we did not find any commercial simulator fulfilling the aforementioned conditions with such a reduced sampling time that was open-source and internally customized. Next, the architecture of the proposed simulator is presented, which would allow testing chlorination control algorithms for treatment plants (DWTP), fulfilling the following features:


One of the main advantages of working with a model implemented in an open source language is the possibility to start with a simple approximation of the process and then gradually refine (programming) the model as the understanding of the process improves. The continuous process of refining leads to achieving a good (accurate enough) approximation to solve the problem of the chlorination tank simulator fulfilling the previously mentioned requirements.

The contribution of this paper is the design and implementation of such a simulator. The architecture depicted in Figure 2 is the key of the contribution, and the starting point for the implementation. With the proposal of this architecture, we will fulfil the four requirements described above. The simulator integrates the equations to describe a hydraulic model and a transport model of diluted species—in this case, the chorine reaction.

**Figure 2.** General architecture of the simulation system for tank chlorination.

The tank behavior (*Simulator Tank* block) was implemented by means of a Simulink S-function, which works in real time, communicating the simulator with the control algorithm implemented in the PLC. The communication of the control loop is the OPC platform. This communication will allow for the importing of SCADA real data from the plant to compare and validate the simulator behavior. In addition to the control variable, the output of the controller transfers the water speed and the level of the tank to the simulator. Conversely, the simulator provides two outputs: (1) the response of the simulated plant to the input of the control algorithm, and (2) the error between the measured chlorine dose and the necessary dose. On the other hand, the output of the controller provides a set of variables such as temperature, speed, or the distribution of flows between osmosis and carbon filters, among others, which allow the chlorine decay to be estimated using the *Decay Function* block. The *Comsol CFD (3D)* block provides a table of diffusion coefficients for locations at different distances of the chlorine output meter of the treated tank in an offline process. In the next sections, we present how this design was implemented and validated, showing relevant results compared with real data.

## *2.2. Description of Sant Joan Despí DWTP*

Since its inception in 1955, the DWTP of Sant Joan Despí (Barcelona, Spain) has aimed to purify the surface waters captured from the Llobregat River and groundwater from the delta aquifer of the river, today supplying drinking water to more than half of the inhabitants (1,900,000) of the city of Barcelona, as well as the southwest conurbation.

The DWTP in Sant Joan Despí incorporates the latest technology in drinking water treatment. The schematic of the plant is shown in Figure 3. Since the water is captured in the Llobregat river, the following processes are applied: surface water capture and desanding; peroxidation; pumping water to the decantation tanks (tag 5 in Figure 3); filtration by sand; pumping by means of four Archimedes screws, wherein the water is raised so that it can follow the process by its own gravity; ozonation; filtration by activated carbon filters; and reverse osmosis, to finally arrive at the mixing and chlorination chambers (tags 13, 14 in Figure 3). In these chambers, the water from the two treatment lines is mixed and chlorination takes place, which ensures the removal of almost all the ammonium content remaining in the water. The study of this paper concentrated mainly on these tanks. Finally, two pumping stations inject water to the distribution network at different levels to be supplied.

**Figure 3.** Drinking water treatment plant (DWTP) at Sant Joan Despí schematic.

The DWTP has a series of processes prior to chlorination in the effluent that favor, to a greater or lesser extent, the elimination of organic matter of natural source and those derived from human activities.

One of the last processes of the plant is chlorination, before driving the treated water to the tanks distributed throughout the city of Barcelona. In [11], there is a good explanation of how the Barcelona water network is implemented, not only physically, but also with an extensive set of sensors to obtain chlorine concentration at the storage tanks distributed along the city, with a total amount of 200 sensors for a distribution network of 4600 km. These tanks, close to the final destination of drinking water (domestic homes), have sensors that measure the lack of chlorine during its transportation from the DWTP to the storage tanks, and a last rechlorination is done to ensure a concentration in the range of 0.8 to 1 ppm when reaching the taps.

To better understand the power of the DWTP, Table 1 contains a set of parameters of water taken at different points in the plant, specifically, the sampling points located: (i) at the entrance of the plant, the intake; (ii) at the entrance of the mixing chamber in the chlorination tanks; and (iii) at the output of chlorination tanks when the drinking water is ready to be transported and distributed through the water network. The chlorination process was performed appropriately because the ammonium concentration dropped to 0 ppm. This is one of the most critical parameters to control to avoid, and the rest of parameters were expected to have those values after all the processes involved in the plant.

In the following sections, the chlorine disinfection strategy in both tanks of the plant is described. Before the process was carried out automatically by the control and instrumentation electronics and the supervision system, we performed the process manually. The facilities and equipment that compose the chlorine regulation system are detailed as well. It should be borne in mind, however, that the validation step described in this article was carried out in the first tank (tank 1) where the ammonium

appears and can be detected. The same process would be totally valid when applied to the second tank (tank 2).


**Table 1.** Parameters of water at different sampling points in the DWTP.

\* Chlorination tank inlet; \*\* to distribution network; \*\*\* this value is due to dosage at the mixing chamber.

#### *2.3. Disinfection Strategy*

Under normal operation, the water previously treated in the different processes of organic matter elimination enters the disinfection phase by combining the diluted chlorine gas and the water to be treated in two large tanks. The retention time required for disinfection is a delay time that makes feedback control a difficult task. The disinfection strategy is based on covering the basic demand of chlorine in the first tank, that is, the needs for chlorine demand are met as a result of the consumption of certain existing organic and inorganic matter. Then, in the second tank, the necessary chlorine is added, that is, free chlorine, to cover losses (by chemical reactions, not leakages) by its transportation. In the same way, the added chlorine needs to reach a certain level to comply with the legal regulations in the effluent (see references above for regulations), with a maximum value of 1.0 mg/L for a good taste and a minimum of 0.2 mg/L to ensure the death of pathogens, bacteria, and viruses. Currently the concentration values required at the exit of the DWTP are usually between 0.8 and 1.2 mg/L of free chlorine.

As shown in Figure 4, the dosing stages to be carried out in the first tank (tank 1) would be between zones 1 and 3. The water to be treated would come into contact with the chlorine in the first tank (zone 1) to react quickly with inorganic matter such as Fe++, Mn++, H2S, or organic matter. Once the initial demand for chlorine has been satisfied as dosing continues, the creation of monochloroamines, dichloroamines, and trichloroamines (zone 2) begins due to the combination of free chlorine with ammonium or its derivatives. In drinking water, monochloroamines are deliberately formed by the reaction of chlorine in an aqueous solution with the added ammonium ion or non-dissociated ammonia in water [12]. In zone 2, total chlorine is a combination of free chlorine and its different organic compounds, and its disinfection effectiveness is limited when total chlorine is only formed by free chlorine. If the dosage continues increasing (zone 3) when all the nitrogen-derived species have oxidized, the reaction reaches a breakpoint of special importance in the chlorination process because beyond this breakpoint the entire dosage will be fully effective in terms of disinfection.

From the beginning of zone 4, all the dosed chlorine becomes residual free chlorine, becoming the ideal point of entry to the tank 2. Specifically, the disinfection strategy that is followed in the DWTP is to raise the outlet from tank 1 to 0.2 ppm to ensure that tank 2 operates only with free chlorine.

**Figure 4.** Transition in the dosing areas of the chlorination in tank 1 and inlet of tank 2.

#### *2.4. Tanks and Equipment*

As previously mentioned, the disinfection system has two contact tanks, tank 1 and 2, with a capacity of 10,000 m<sup>3</sup> each. These tanks are hydraulically coupled and, in case of an anomaly in any of them, the DWTP can treat water only with the available tank (only for a limited time of operation). The water flow to be treated in tank 1 (Figure 5) arrives at a small chamber where water treated with a traditional treatment by granulated active carbon filters (GACF) is mixed with water treated with a reverse osmosis (RO) treatment. At the exit of both treatment processes, there are flow meters that allow sensing the flow coming from both the osmosis and the carbon filters. The sum of these two flows will be the variable of total flow to be treated, as well as the input for the control system and the simulator presented in this article. The actual flow range treated by the plant is between 5.5 m<sup>3</sup> /s at maximum performance and a minimum flow in certain operating situations of 1.0 m<sup>3</sup> /s. In normal regime, the distribution by GACF and by RO is 50%, however, this distribution may change in a few hours due to the operating conditions of the DWTP. The theoretical residence time of water in the one tank can vary between 23 min and 55 s for a flow rate of 5.5 m<sup>3</sup> /s up to 2 h and 11 min for a flow rate of 1.0 m<sup>3</sup> /s.

**Figure 5.** Plant view of chlorination tank 1 (measured in meters).

In the inlet of tank 1, the mixing chamber, we installed an ammonium meter (An1), allowing the control system to acquire the amount of ammonium contained in the water to be treated, with the first chlorinated water dosage point (1) being located a few meters apart downstream of the entrance.

It is essential for the control system to know the delay times between the production of chlorinated water and the dosage in point 1 since they are dead times that must be compensated and considered in the calculation of the dose. At a distance of 87 m from the dosing point 1, there is the first residual chlorine meter (2). This distance between the dosing point and the measurement of free chlorine is considered sufficient in most cases for the chlorine to react with organic matter and other compounds. Depending on the velocity of the water and the amount of ammonium, sometimes the measurement is no longer representative for errors in the meter and for this reason there is a second meter (3) located at the exit of the tank. As shown in Figure 5, the tank is divided into two parts and forces the treated water to travel twice the distance so that the reaction is complete and the contact time is sufficient for an effective disinfection. In the center of the tank, there is an ultrasonic level meter with a height range of 0–4 m.

The analytical instrumentation provided by the control system is composed of a colorimetric ammonium-type analyzer at the entrance of the tank and two amperometric meters that measure the residual chlorine at the entrance and exit of the two tanks. The response of the amiometric residual chlorine meters for combined chlorine samples has been quite satisfactory once the reaction time of the ammonium with the free chlorine is guaranteed. To measure the flow at the output of each of the 20 GACF and the 10 RO frames, we used electromagnetic flow meters. It is important to parameterize these flowmeters at the level of integration of the measure to dampen possible sudden changes in the reading of the control system. Finally, an ultrasonic level meter is available in both tanks that allow us to know the variations of the volume stored in each tank on the basis of the knowledge of their dimensions.

#### **3. Mathematical Process Modeling**

In this section, the graybox model implemented is presented. This model allows for a simulation of the chlorination process that describes the behavior of the fluid by the advection and diffusion equations, as well as a chlorine-demand estimation function to model the reaction.

The chlorine dosing system in contact tanks 1 and 2 can be characterized as a distributed parameter system that takes place in both coupled reactors. The chlorine dosed at the beginning of each tank (reactor) is subject to a first-order decay which allows the mass balance shown in Equation (1). This partial differential equation (PDE) represents the model of a solute in a turbulent flow, in 3D:

$$\frac{\partial c}{\partial t} = D\_{t-x} \frac{\partial^2 c}{\partial x^2} + D\_{t-y} \frac{\partial^2 c}{\partial y^2} + D\_{t-z} \frac{\partial^2 c}{\partial z^2} - \mathcal{U} \frac{\partial c}{\partial x} - V \frac{\partial c}{\partial y} - \mathcal{W} \frac{\partial c}{\partial z} \tag{1}$$

where *c* is the chlorine concentration; *D t-x,y,z* is the dispersion; and *U*, *V*, and *W* are the water velocities in the coordinate axes *x*, *y*, and *z*, respectively. Simplifying it to a 1D equation and adding the reactive term, we obtained the following equation:

$$\frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial \mathbf{x}^2} - \mathbf{U} \frac{\partial c}{\partial \mathbf{x}} - kc \tag{2}$$

where *U* is the water velocity in the tank (m/s). Therefore, Equation (2) is the advection–diffusion–reaction equation (ADR) and belongs to the group of partial differential equations of parabolic type.

The simulator developed and presented in this article is based on the discretization of Equation (2) [13], treating the chlorine decay part (*kc*). The diffusion coefficient values were obtained from the Comsol CFD Software on the basis of the tank and flow characteristics, treated as a RANS *k*-ε turbulence model.

The proposed simulator, despite working in 1D, collects all the necessary information for the characterization of the process and allows the control system to be tested as if it were actually being dosed in the tank. In the following sections, we detail the more relevant aspects and considerations that were considered at the time of the design, as well as the implementation of the simulator related to advection transport, flow characteristics and the parameterization of the chlorine reaction, and diffusion of the chlorine through the tank.

In the following sections, the turbulent flow regime (Reynolds number) is verified, and the Péclet number is estimated in order to assess the influence of the convective versus diffusive effect. Finally, the Schmidt number chosen for this turbulent flow model is presented.

#### *3.1. Flow Characteristics*

The flow in the tank under study corresponds to a three-dimensional speed distribution. This is a non-stationary flow, considering the fluctuation of the flow in the tank and not its uniformity due to its geometry considering the characteristic velocity profile of the flow in open channels.

The consideration of constant velocity throughout the domain would be a relatively correct simplification if the chlorine were evenly distributed at all times. However, in principle, this hypothesis is not acceptable due to the characteristic velocity profile for this type of geometry. In addition, the geometry of the studied domain should be considered where changes occur in the direction of flow. This situation even drives away the actual behavior of the fluid from the idealized one-dimensional flow situation.

On the other hand, we faced a turbulent flow, according to the Reynolds number range, Equation (3), for the featured width of the tank L(width) = 8.2 m and according to the properties of the water at normal working temperature, with a velocity for maximum and minimum operations points [0.052 . . . 0.128]m/s being

$$Re = \frac{\rho vL}{\mu} = \frac{1000 \times v \times 8.2}{0.001} \tag{3}$$

where ρ is the density of the fluid (gr/L), *v* the velocity, *L* the tank width, and µ the dynamic viscosity of the fluid. The water velocity range was obtained on the basis of the flow and the elapsed time to cross the tank of known dimensions. The minimum and maximum values for the Reynolds number were 4.29 <sup>×</sup> <sup>10</sup><sup>5</sup> and 1.06 <sup>×</sup> <sup>10</sup><sup>6</sup> , respectively, defining a range of values within the turbulent regime. In [14], a comparative study is presented between different configurations of contact tanks for a disinfection process and the calculation of the Reynolds number, using as a parameter, in this case, the kinematic viscosity of the water as opposed to the dynamic viscosity such as in this research.

#### *3.2. Di*ff*usion Coe*ffi*cient*

The diffusion coefficient depends largely on the nature of the particles, the solvent, the temperature, and the viscosity of the solvent. In this case, neither the type of particles nor the solvent is modified during the process. On the other hand, temperature variations are not important with respect to diffusion and, as a consequence, there are also no significant variations in the viscosity of the solution. Considering all these conditions, we can affirm that the diffusion coefficient will remain constant.

In general, if there is an active flow, as presented in this paper, the diffusion effect is negligible. However, to determine the relationship between the advective and diffusive terms, the Péclet number *Pe*, Equation (4), was estimated:

$$Pe = \frac{L \times v}{D} \tag{4}$$

where *L* is the distance between concentration measurement points, *v* the fluid velocity, and *D* (m/s 2 ) the diffusion coefficient.

In this study case, the Comsol CFD software was used to obtain the value of the diffusion coefficient of the particles from the characterization of the tank and different flow rates giving a range of 0.0055 to 0.013 m<sup>2</sup> /s for water flow in the range between 1.5 m<sup>3</sup> /s and 3.7 m<sup>3</sup> /s, respectively. Similar to this strategy, in [15], CFD software was used to obtain the dissemination number of the treated tank.

For the calculation of the minimum and maximum values of the number of Péclet, according to the extreme values of the velocity, we took the most unfavorable value of the diffusion coefficient (*D* = 0.013 m<sup>2</sup> /s). Considering these conditions, the Péclet number oscillated between 71.92 and 1831.0, both higher than 1. These values suggest that the diffusive term was negligible compared to the effects

of advection, and since the diffusion can be neglected for *Pe* >> 1 under these conditions, the effects of advection exceeded those of diffusion in determining the overall mass flow.

#### *3.3. Turbulence Model*

The turbulence model chosen for the CFD study is the R-based k-epsilon (*k*-ε) model, which is a two-equation model that provides a general description of the turbulence through turbulent kinetic energy (*k*) and dissipation (ε). The *k*-ε turbulence model is widely used to model flow behavior in chlorination tanks of these characteristics [16].

The Schmidt number (*Sc*), Equation (5), is a relevant parameter in the configuration of the model used in this study, which is used in RANS *k*-ε to avoid the resolution of the boundary layer. The *k*-ε model predicts turbulent viscosity thanks to the Schmidt number:

$$S\_{\mathcal{E}} = \frac{\mu}{\rho D\_{\mathcal{I}}} \tag{5}$$

where µ is the dynamic viscosity of the fluid, ρ the density of the fluid, and *D<sup>t</sup>* is mass diffusivity of the fluid—a value that depends on the fluid, water plus chlorine in this case.

The Schmidt number is an empirical constant with typical values between 0.1 and 1; in the present study giving a value of 0.7, following the criteria of many other documented works on water treatment and contact tanks of similar characteristics [17–20]. This parameter is relatively insensitive to the properties of the molecular fluid (the particular value obtained from other experiments can be used despite the differences between the simulation domain). It represents a significant parameter for fully developed turbulent flows, and it is considered that, for the present case, the average Reynolds number is in the transition to the turbulence zone with moderate values, and thus it was not a key parameter for the present case. A sensitivity analysis of the *Sc* was carried out in the present study, showing the little relevance of the changes in its value close to 0.7; between 0.5 and 0.9 the concentration variations obtained in the effluent were around ±0.01 ppm. Finally, with respect to *Sc*, it should be borne in mind that the experimental determination of this parameter was not within the scope of this project.

#### *3.4. Chlorine Reaction through the Tank*

The expression in Equation (2) is a generalization that must be developed with more precision for the present case. In [21,22], different models are presented, describing a chlorine decay more adjusted to the reality of the studied plant. Considering the characteristics of the water that reaches the exit of the DWTP, a significant part of the chlorine (fraction *f* in Equation (6)) reacts quickly with the existing matter (organic and inorganic) and the rest (fraction 1−*f*) of the matter reacts with the remaining chlorine fraction. Therefore, a combined first-order model according to [21], plus a combination of first and second order presented in [22], yields Equation (6):

$$\frac{\partial \mathcal{C}}{\partial t} = D \frac{\partial^2 \mathcal{C}}{\partial \mathbf{x}^2} - \mathcal{U} \frac{\partial \mathcal{C}}{\partial \mathbf{x}} - k\_{\mathcal{R}} \mathbf{c}(f) - k\_r \mathbf{c}^2 (1 - f) - k\_{\mathcal{S}} \mathbf{c} (1 - f) \tag{6}$$

where *k<sup>R</sup>* is the decay coefficient for the rapid reaction, and *k<sup>r</sup>* and *k<sup>s</sup>* would be the decay constants of rapid and slow reaction for the remaining chlorine fraction with different characteristics to those that react following *kR*. Let *f* be the fraction of chlorine that reacts quickly in response to the decay constant *kR*.

It has been proven by laboratory analysis that, for treated waters without ammonia, the demand for chlorine varies between 0.2 ppm and 0.5 ppm depending on whether its source is from the traditional treatment of granular active carbon filters (GACF) or osmosis (RO). These reactions take place between 5 to 10 min after the chlorine is exposed to the water and follows a decay according to the constant *kR*. The demand for chlorine is very high, and *k<sup>R</sup>* varies significantly when the existence of ammonium occurs, with values between 7 and 10 times the amount of NH4-N detected by the ammonium

meter (An1). Constants *<sup>k</sup><sup>r</sup>* and *<sup>k</sup><sup>s</sup>* were estimated in [23], with values between 3.3934 <sup>×</sup> <sup>10</sup>−<sup>5</sup> (s−<sup>1</sup> ) and 5.4259 <sup>×</sup> <sup>10</sup>−<sup>7</sup> (s−<sup>1</sup> ), respectively. These values (*k<sup>r</sup>* and *ks*) are significant because of the importance in the parameterization of the proposed simulator. With those values of *k<sup>r</sup>* and *k<sup>s</sup>* , and knowing that the demand for chlorine in its rapid reaction stage is between 0.2 ppm and 1.0 ppm and that the fraction that reacts with the chlorine is around 40%, then *<sup>k</sup><sup>R</sup>* is between 3.046 <sup>×</sup> <sup>10</sup>−<sup>4</sup> (s−<sup>1</sup> ) y 14.05 <sup>×</sup> <sup>10</sup>−<sup>4</sup> (s−<sup>1</sup> ), according to experimental data.

Figure 6 shows the simulator response to a step in the dosage of 1 ppm for a flow of 1.5 m<sup>3</sup> /s parameterized with *k<sup>r</sup>* and *k<sup>s</sup>* , and different values of *k<sup>R</sup>* with a reaction fraction of *f* = 0.4 (40%).

**Figure 6.** Step response for a dose of 1 ppm of chlorine for different values of *k<sup>R</sup>* (0.0001, 0.0005, 0.0009, 0.0013 s−<sup>1</sup> ).

ε − − Figure 7 shows a view of one of the studies carried out in steady state with the Comsol CFD on the diffusion of chlorine along the contact tank. Specifically, on the basis of a color scale, the diffusion of chlorine along the tank is represented for a turbulent flow rate RANS *k*-ε with a flow rate of 2.7 m<sup>3</sup> /s, with a decay constant of 9.0 <sup>×</sup> <sup>10</sup>−4 (s−<sup>1</sup> ), a Schmidt number of 0.7, and an initial dosing step of 2 ppm of chlorine. ε − −

− − **Figure 7.** Comsol simulation for decay constant of *<sup>k</sup><sup>R</sup>* <sup>=</sup> 9.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> (s−<sup>1</sup> ).

− −

#### **4. Discretization**

0= ଶ <sup>ଶ</sup> −

0= ଶ <sup>ଶ</sup> −

−

−

The partial derivative equation (PDE) presented in Equation (6) in continuous space was reduced to a second-order ordinary differential equation (ODE), Equation (7), in the discretization of the

ଶ

ଶ

(1 − ) − ௦(1 − )

(1 − ) − ௦(1 − )

− ோ() −

− ோ() −

simulator. It was considered that the type of simulator flow belongs to the "Plug Flow" pattern where the fluid circulates through the tank evenly and along parallel paths from the entrance to the exit of the tank.

$$0 = D\frac{d^2\mathcal{C}}{dx^2} - \mathcal{U}\frac{d\mathcal{C}}{dx} - k\_{\mathcal{R}}c(f) - k\_{\mathcal{r}}c^2(1-f) - k\_{\mathcal{s}}c(1-f) \tag{7}$$

Among the different methods of discretization for the type of equations that govern the behavior of the simulated system process (Quick, Upwind, etc.), we chose a central finite differences equation, which was adapted in an acceptable way for the values of flow, diffusive characteristics, and Péclet numbers of the present problem [24]. In [25], different discretization methods were presented and discussed to address advection, diffusion, and reaction problems in contact tanks; in [26], a discretization scheme was presented for storage tanks in transport processes and distribution of the same type of water treated by the DWTP of this article.

When applying central finite differences by replacing the first and second derivative in Equation (7), then Equation (8) is obtained:

$$0 = D \frac{\mathbb{C}\_{i+1} - 2\mathbb{C}\_i + \mathbb{C}\_{i-1}}{\mathfrak{x}^2} - \mathbb{U} \frac{\mathbb{C}\_{i+1} - \mathbb{C}\_{i-1}}{2\mathfrak{x}} - k\_{\mathbb{R}} \mathfrak{f}(f) - k\_r c\_i^2 (1 - \mathfrak{f}) - k\_s c\_i (1 - \mathfrak{f}) \tag{8}$$

for a diffusion coefficient D and a chlorine concentration c at different instants of time.

Figure 8 shows a diagram of the diffusion effect in chlorine concentration. It supposes a tank of length *L* with constant increments of space on the mesh (∆x) at different periods of time. An increment of <sup>∆</sup>x = *L*/(*n* <sup>−</sup> 1) was considered, with *n* being the number of mesh divisions. To avoid the effects of dynamic instability, the increment values of x fulfilling *x* ≤ 2*D <sup>v</sup>* were restricted, and the stability criterium *t* ≤ (*x*) 2 2*D*+*kRx* <sup>2</sup> was applied according to [27]. Applying the stability restriction on the spacing <sup>∆</sup>x in the most unfavourable case, a value of <sup>∆</sup><sup>x</sup> <sup>≤</sup> 0.2 m was obtained and <sup>∆</sup><sup>t</sup> <sup>≤</sup> 2.16 s. In the simulation, knowing that the actual dosage analyzer 2 (An2) is located at 83 m from the tank 1 entrance, a value ∆x = 0.127 m and ∆t = 1 s with a mesh of *n* = 685 slots was used, in accordance with a trade-off to minimize the error with a computing time that would allow the simulation in real time. Δ ∆ − ∆ ≤ ଶ ௩ ∆ ≤ (∆௫) మ ଶାೃ∆௫<sup>మ</sup> Δ Δ ≤ Δ ≤ Δ Δ

**Figure 8.** Diagram representing the diffusion effect respect time and space (not at scale) in the tank for a water inflow *Fin*. At different instant times (t<sup>0</sup> , ti , tj . . . ), the concentration of chlorine is different in the tank (c<sup>0</sup> , c<sup>i</sup> , cj . . . ) depending on the distance to the tank inlet (function of *x*).

#### **5. Simulator Implementation**

For the implementation of the simulator, the Matlab/Simulink tool was used using the C language to encode the functions that describe the behavior of the model. The simulation code was programmed and encapsulated in an S-Function (Level 2) in Simulink environment and the interface was adapted to be able to contrast data from the output of the supervision and control software (SCADA) with the output of the simulation. The proposed simulator could be configured for any programmable logic controller (PLC) and SCADA controller with an OPC interface (OLE for process control, object linking, and embedding). One of the advantages of using a S-Function is that it allows for the creation of a general-purpose block that can be used in all the iterations of the model, varying the parameters with each instance of the block and thus interacting in real time with the controller (PLC). The simulator implements the central finite differences for the discretized model in Equation (8) [28].

Simulation main screen (in Simulink environment) for the interaction scheme between the simulator and the PLC controller according to the simulator architecture (Figure 2) is shown in Figure 9. The data from the controller are evaluated by the *estimateCd* function that calculates the chlorine demand and, thus, the value of the rapid decay coefficient (*kR*) for each iteration. Depending on the velocity of the flow, we determined the diffusion coefficients that are part of the input variables of the simulator, together with the dose, *kR*, and velocity. The designed interface allows for the setting of the distance to the reading point, the rapid reaction fraction for *kR*, and constant values *k<sup>r</sup>* and *k<sup>s</sup>* .

**Figure 9.** Connection between the simulator and a programmable logic controller (PLC) in the OPC (object linking and embedding (OLE) for process control) platform.

The simulator comprises two working zones. If there is less than 0.02 ppm of ammonium in the treated water (sensed with analyzer An1), the chlorine demand is calculated on the basis of the temperature and the source of treated water. We developed a classification of the chlorine demand on the basis of the water temperature (Temp) and the percentage of water passing through the carbon filter treatment (TPGacf) with respect to the total treated flow. In the tree, each child node represents a chlorine demand value based on the classification by temperature and source and percentage of treated water. This study was generated through taking real data from the last 5 years in the plant. Figure 10 (left) shows the decision tree for classified data. The study was implemented with RStudio. For values greater than 0.02 ppm in the ammonium concentration, the meter yields accurate and reliable results, having proven experimentally that there is a linear correlation (greater than 0.93) between the ammonium concentration and the chlorine demand for treated water (see Figure 10, right).

The validation of the simulator is an important part of this work, because it allows for the measurement of the goodness of such a simulator in its use in forecasting different situations in the real plant in front unknown events that can happen in the future. Results are shown in the next section, but here the implementation of how the simulation is validated is presented.

**Figure 10. Left**: decision tree for chlorine demand; **right**: correlation between ammonium concentration and chlorine demand (about *R* = 0.93).

The procedure to compare real data obtained with the SCADA and simulated data is shown in Figure 11. Two blocks (*Tank\_secA* and *Tank\_secB*) are instantiated in the test, one for each installed analyzer, An2 and An3, at 87 and 163 m away from the dosing point, respectively. The level of the tank (*Level*) and the flow rate (*TotalFlow*) are used to calculate the velocities in sections A and B of the tank. The function blocks (*An2 Sample* and *An3 Sample*) simulate the measurement behavior of chlorine analyzers with periodic samples every 5 min. The *Chlorine Dosage System* function block simulates the process of injecting chlorine gas into water from its generation point and its delay in transport to the contact tank, in order to be more realistic in the simulation procedure.

**Figure 11.** Simulator block diagram to compare simulation data and real data from SCADA.

At present, the chlorination stage is performed in a semi-automatic regime. From the plant control room, the operators check the level of chlorine concentration on the basis of the sensors' readings, and they order PLC (programmable logic controllers) to open the valves of chlorine bottles, which adds a specific dosage to reach the demand at the control point (located as indicated in Figure 5).

#### **6. Results: Validation of the Simulator**

The validation of the simulator was divided into two stages with specific objectives. In a first stage, data from the Comsol CFD, considered as a reference computational fluid dynamics software, were contrasted with respect to the output of the simulator for conditions of advective flow, diffusion, and decay of chlorine. The objective of this phase was to verify the correct implementation and validation of the advection–diffusion–reaction model in the Simulink S-Function. In a second stage, the *estimateCD* function was parameterized and implemented for the characteristics of the water treated by the DWTP, and data from the SCADA database were validated on the basis of this parameterization.

Table 2 shows a comparison of values obtained from Comsol CFD with respect to those provided by the simulator for different decay constants and a flow rate of 3.7 m<sup>3</sup> /s. In any case, the relative error in the sample is larger than 2% considering that chlorine meters have an absolute error of 0.04ppm. This fact implies that the error made by the simulator would widely meet the requirements in terms of accuracy.


**Table 2.** Comparison of chlorine decay between the simulator and the Comsol computational fluid dynamics (CFD). −

Figure 12 shows graphically a comparison of data obtained with Comsol CFD and the simulator, with different water flows between 1.50 and 3.70 m<sup>3</sup> /s. On the basis of these results, the simulator gave excellent results.

**Figure 12.** Results comparison between CFD Comsol data and simulator data for different flow rates (1.50, 2.00, 2.70, and 3.70 m<sup>3</sup> /s).

However, the most important validation of simulator was carried out with real data; this is perhaps the most effective way to demonstrate the power of a simulator. In this case, the experiments were performed with real data obtained from the plant in real situations. Most of the time, the plant has neither variation on the input flow nor the effluent where the water comes from, and thus the chlorine dosage is constant. Such situations are easy to model. However, when weather conditions are varying, such as in episodes of flooding or big storms, or even in large episodes of draught when the river flow in very low, the situation in the plant becomes exceptional. In such periods, the chlorine dosage has to abruptly change due the bad conditions of input flow to the plant. We focused our attention upon such

days, and obviously to the associated data recorder during those extreme episodes. Together with the flow, the dosage is already recorded as in a log file. Therefore, it is easy to reproduce the conditions that generated a specific behavior in the plant, and those conditions were used as input in the simulator to recreate the real situation that occurred in the plant. The output of the simulator was compared with the chlorine concentration in the measurement point (output of the plant) as a consequence of the dosage that was injected in the inflow water. This comparison study was carried out in all the operating points of the plant, that is, at the different input flows of 1.5, 2.0, 2.7, and 3.7 m<sup>3</sup> /s.

As can be seen in Figure 13, the plant behavior refers to the chlorine concentration at the end of the tank, that is, the concentration of chlorine in parts per million ready to enter the pipes for distribution. The model is the execution result of the simulator as the chlorine concentration at the same point of the plant, and therefore plots in Figure 13 compare both concentration in the same location of the tank. As it can be observed, the simulator behavior worked exceptionally well following the reality in front of different input flow situations. This qualitative assessment aside, the variations between real data and simulated data are quantified in Figure 14.

**Commented [M15]:** Please change the comma to a **Figure 13.** Data comparison between real data from SCADA and simulator data for different flow rates (1.5, 2.0, 2.7, and 3.7 m<sup>3</sup> /s).

decimal point in the figure 13 **Commented [UdW16R15]:** Done

**Figure 14.** Plot of residuals for different flow rate operating points (1.5, 2.0, 2.7, and 3.7 m<sup>3</sup> /s).

Defining a residual value (or residuals with sign) as the difference between the observed values and the values predicted by the model (estimated values), Figure 14 shows a graph of the residual values obtained in the different operating points. For all flow values, a distribution of bell-shaped residuals can be seen, except for flows close to 1.5 m<sup>3</sup> /s, which is multimodal. At large flow rates with high velocity (above from 1.5 m<sup>3</sup> /s), it is easier to predict errors and residuals closer to zero and with less dispersion, because the dispersion, either multimodal or unimodal, increases when the flow drops. One explanation is that the simulation treats the problem as a case of one-dimensional flow (and Plug Flow), and it does not consider that the velocity is not uniform throughout the domain. Not all the dosed chlorine moves at the same velocity in the tanks and it is also possible that there are fluid recirculations in the studied domain. As a consequence, at the measurement point, the concentration is detected with a time lag with respect to the simulation, and at low flow rates this fact is greatly accentuated. Finally, the chlorine analyzer error of ± 5% over 5.00 ppm must be considered and, therefore, those simulation errors can be perfectly accepted. The abscissa axis indicates the error in concentration of chlorine between real and simulated data. The ordinate axis represents the frequency of such residuals appears, that is, a kind of histogram of number of appearances. Note that the frequency is in thousands because we studied large amounts of data, searching for exceptional situations in the plant, and then for every particular study the sequence is for 241 min, sampled at 1 s, the sampling time of the designed simulator.

As it can be seen in Table 3, the mean square error (MSE) are collected as a measure of the quality of the simulator. In this case, the difference between the real data and the simulated data is measured as:

$$MSE = \frac{1}{N} \sum\_{t=1}^{N} (c\_{\text{real\\_data}} - c\_{\text{simulated\\_data}})^2 \tag{9}$$

that is, the sum of the square differences between the chlorine concentration in the real plant and the result of the simulator in front the same dosage and flow. It is worthwhile to note that variable *t*, the time, has a step of 1 s.

**Table 3.** Mean square error (MSE) between the observed values (real data) and the estimated values (simulated data).


#### **7. Conclusions**

In this article, we proposed the design and implementation of a simulator that allows for the verification of control methodologies in a simple and visual way before its implementation directly in a real drinking water treatment plant DWTP. The proposed simulator is specific for processes that involve a hydraulic model where transport of diluted species can be simulated. The simulator provides simplicity, easy connection to plant control equipment (with OPC platforms), and reliability. The connection between the specific hardware of the plant and the simulator operated in a satisfactory manner, allowing a data interchange at a sample rate of 1 s, time enough in this kind of phenomenon. The simulator was validated with a more complex simulator (Comsol CFD) unable to operate with the plant equipment. The result of such a comparison was very satisfactory, giving an error less than 2% in the worst case. Moreover, and a crucial test, the simulator was validated with real data acquired under extreme circumstances in tough periods in terms of water effluents arriving to the DWTP. The results were completely satisfactory, with a good performance. Those good results confirmed that the reduction of the model from partial derivatives equations to ordinary differential equations was correctly performed, and thus the model is reliable enough to be confident in its functioning in the plant. After these results, we are confident in the use of the simulator as a forecast tool in a real plant. **Author Contributions:** Conceptualization, J.G. and A.G.; methodology, J.G. and Y.B.; software, J.G. and H.M.; validation, H.M. and Y.B.; formal analysis, A.G. and H.M.; investigation, J.G. and Y.B.; resources, A.G. and H.M.; data curation, J.G. and Y.B.; writing—original draft preparation, J.G.; writing—review and editing, A.G.; visualization, H.M. and Y.B.; supervision, A.G. All authors have read and agreed to the published version of the manuscript.

**Acknowledgments:** We would like to thank Francisco Luque Montilla; David Ibarra from Aquated, Inc.; Aigües de Barcelona analytic laboratory staff; Ricardo Torres, UPC; and Mercedes Garcia, Education Consortium Barcelona.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
