**Energy Saving at Cities**

Printed Edition of the Special Issue Published in *Energies* Francisco Manzano Agugliaro and Alberto Jesús Perea Moreno Edited by

www.mdpi.com/journal/energies

**Energy Saving at Cities**

## **Energy Saving at Cities**

Editors

**Francisco Manzano Agugliaro Alberto Jes ´us Perea Moreno**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Francisco Manzano Agugliaro University of Almeria Spain

Alberto Jesus´ Perea Moreno University of Cordoba Spain

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Energies* (ISSN 1996-1073) (available at: https://www.mdpi.com/journal/energies/special issues/ energy saving cities).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-760-3 (Hbk) ISBN 978-3-03936-761-0 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editors**

**Francisco Manzano Agugliaro**, full professor at the Engineering Department in the University of Almeria (Spain), received his M.S. as Agricultural Engineer and Ph.D. in Geomatics at the University of Cordoba (Spain). He has published over 160 papers in JCR journals (https://orcid.org/0000-0002-0085-030X), H-index 26. His main interests are energy, sustainability, scientometrics, water, and engineering. He has served as a Supervisor of 25 Ph.D. Theses; Vice Dean of Engineering Faculty (2001–2004); Director of Central Research Services (2016–2019); Ph.D. Program Coordinator for Environmental Engineering (2000–2012), Greenhouse Technology, Industrial and Environmental Engineering (from 2010); and General Manager of Infrastructures (from 2019) at University of Almeria. Awards: Top reviewers in Cross-Field—September 2019 (Web of Science), 2019 Outstanding Reviewer Award (Energies), 2019 *Sustainability* Best Paper Awards.

**Alberto Jes ´us Perea Moreno**, associate professor at the Department of Applied Physics, Radiology, and Physical Medicine in the University of Cordoba (Spain), received his M.S. as Agricultural Engineer and Ph.D. in Geomatics at the University of Cordoba (Spain). He has published over 40 papers in JCR journals (http://orcid.org/0000-0002-3196-7033), H-index 12. His main interests are renewable energy, energy saving, biomass, sustainability and remote sensing. He served as the secretary of the Applied Physics Department (2017–2019) at University of Cordoba. Awards: 2019 *Sustainability* Best Paper Awards.

## **Preface to "Energy Saving at Cities"**

The use of renewable energies and energy saving, and efficiency are the needs of a global society. According to the latest estimates, global energy demand could triple by 2050, and by then, 70% of the world's population will live in cities. Cities are currently responsible for 80% of greenhouse gas emissions, so they have a key role to play in shifting toward a sustainable energy future. Cities are threatened not only by overcrowding, but also by rising energy demand, obsolete infrastructure, volatile energy markets, and the effects of climate change.This book aims to advance the contribution of energy saving and the use of renewable energies in order to achieve more sustainable cities.

> **Francisco Manzano Agugliaro, Alberto Jes ´us Perea Moreno** *Editors*

### *Editorial* **Energy Saving at Cities**

#### **Alberto-Jesus Perea-Moreno <sup>1</sup> and Francisco Manzano-Agugliaro 2,\***


Received: 3 July 2020; Accepted: 17 July 2020; Published: 22 July 2020

**Abstract:** The use of renewable energies, energy saving, and efficiency are needs of global society. According to the latest estimates, global energy demand could triple by 2050 and, by then, 70% of the world's population will live in cities. Cities are currently responsible for 80% of greenhouse gas emissions, so they have a key role to play in shifting towards a sustainable energy future. Cities are threatened not only by overcrowding, but also by rising energy demand, obsolete infrastructure, volatile energy markets, and the effects of climate change. This Special Issue aims to advance the contribution of energy saving and the use of renewable energies in order to achieve more sustainable cities.

**Keywords:** energy saving; renewable energy; zero-energy buildings; energy efficiency; sustainability; sustainable transport; PV; energy saving in data processing centers

#### **1. Introduction**

Climate change is increasing due to the anthropogenic emission of greenhouse gases. The majority of these are due to the production and consumption of energy. The challenge for the future cities is the implementation of a mechanism that minimizes the need for injection of new energy resources in them, so that a high level of self-sufficiency should be achieved through the concept of circular economy, thus partially mitigating the impacts of climate change. Using solar energy today is considered to be one of the best solutions that can be installed in buildings to help with this issue.

This Special Issue aims to advance the contribution of energy saving and the use of renewable energies in order to achieve more sustainable cities. This Special Issue sought contributions spanning a broad range of topics related, but not limited to:


#### **2. Publication Statistics**

Details of the call for papers for this Special Issue, regarding the articles submitted being published or rejected, were: 10 articles submitted (100%), 3 articles rejected (33.3%), and 7 articles published (66.6%).

The regional distribution of authors by countries for the published articles is presented in Table 1, in which it is possible to observe 20 authors from six countries. Note that it is usual for an item to be signed by more than one author and for authors to collaborate with others from different affiliations. The mean number of authors per published manuscript was four authors.


**Table 1.** Authors' countries.

#### **3. Authors' A**ffi**liations**

This Special Issue's authors and their first affiliations are reflected in Table 2.


**Table 2.** Authors and affiliations.

#### **4. Topics**

Table 3 summarizes the research carried out by identifying the topics to which they belong, according to the proposed topics in the special issue. It was noted that two "Energy Saving at Cities" topics dominated the rest: "Solar Energy", and "Sustainability".


**Table 3.** Energy Saving at Cities.

**Author Contributions:** The authors all made equal contributions to this article. All authors have read and agreed to the published version of the manuscript.

**Conflicts of Interest:** There is no conflict of interest, the authors state.

#### **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/).

## **Analysis of Research Topics and Scientific Collaborations in Energy Saving Using Bibliometric Techniques and Community Detection**

**Carmen de la Cruz-Lovera 1,\*, Alberto-Jesus Perea-Moreno 1, José Luis de la Cruz-Fernández 1, Francisco G. Montoya 2, Alfredo Alcayde <sup>2</sup> and Francisco Manzano-Agugliaro <sup>2</sup>**


Received: 15 April 2019; Accepted: 24 May 2019; Published: 27 May 2019

**Abstract:** Concern about everything related to energy is increasingly latent in the world and therefore the use of energy saving concepts has been increasing over the past several years. The interest in the subject has allowed a conceptual evolution in the scientific community regarding the understanding of the adequate use of energy. The objective of this work is to determine the contribution made by international institutions to the specialized publications in the area of energy-saving from 1939 to 2018, using Scopus Database API Interface. The methodology followed in this research was to perform a bibliometric analysis of the whole scientific production indexed in Scopus. The world's scientific production has been analysed in the following domains: First the trend over time, types of publications and countries, second, the main subjects and keywords, third, main institutions and their main topics, and fourth, the main journals and proceedings that publish on this topic. Then, these data are presented using community detection algorithms and graph visualization software. With these techniques, it is possible to determine the main areas of research activity as well as to identify the structures of the collaboration network in the field of renewable energy. The results of the work show that the literature in this field have substantially increased during the last 10 years.

**Keywords:** energy-saving; energy conservation; energy utilization; energy efficiency; scientific collaboration

#### **1. Introduction**

The energy problem has become a global issue in a world where science and technology are undergoing a great evolution in a short period of time [1]. Consequently, the development of renewable energy and the use of energy saving concept has been increasing in last years, becoming a research focus around the world [2]. The growing interest in the subject has allowed a conceptual evolution in the scientific community regarding the understanding of the adequate use of energy.

While there is a close relation between the terms "energy saving" and "energy efficiency", there are certain distinctions that need to be clarified. The energetic saving entails a change in the habits of consumption. Sometimes it is enough to eliminate habits that waste energy. On the other hand, energy efficiency is the fact of minimizing the amount of energy needed to satisfy the demand without affecting its quality [3]; it supposes the substitution one machine to another that, with the same benefits, consumes less electricity. Therefore, this does not mean changes in consumption habits. The behaviour of the user remains the same, but less energy is consumed due to the consumption of energy to carry out the same service, which is lower. In Alcott, an extensive review of the literature

on the subject of "sufficiency" can be found [4]. To achieve greater energy sustainability as much as possible, renewable energy and energy saving have been combined [5]. In this way, we will achieve a twofold purpose thanks to a combined strategy of energy saving and energy efficiency, without reducing the quality of life. In fact, it will be maintained, and even increased. These concepts are also making significant advances in the residential construction sectors [6]. Energy saving is not only linked to monetary savings [7], it also contributes to reducing greenhouse gas emissions [8], sustainable development and achieving a sustainable energy supply. However, the reduction of energy demand is not as easy to carry out as is thought and the actions that are imposed are in many cases insufficient to achieve real transformation [9].

Thus, it is necessary to emphasize that, sometimes, many improvements of the energetic efficiency do not reduce the consumption of energy in the amount predicted by the simple engineering models [10]. This is called a rebound effect, since the fact that energy services improve means that they are cheaper and, in turn, that the consumption of these services increases [11]. The Paris Agreement on climate change is an agreement within the framework of the United Nations Framework Convention [12] and energy consumption was the main point to deal with and the creation of an agreement worldwide. To achieve the objectives of this agreement, an outstanding change in the energy sector is necessary, the origin of at least two thirds of greenhouse gas emissions.

The transformation of the electricity sector, with renewable energy, has focused attention on a new debate on the design of the electricity market and electrical safety [13]. However traditional concerns about energy security have not disappeared. In addition, there are also pending issues related to energy access and affordability, global warming and CO2 emissions. Therefore, the massive social concerns can be understood and lead to an increasing demand for ways to improve social and individual energy efficiency.

There remains an urgent demand for ways to improve the energy efficiency of society and people and a need to move increasingly to alternative, low-carbon or non-carbon energy systems. [14]. If we want to achieve sustainable development, energy efficiency becomes a key factor. It is a reality that only a few countries are engaged in research into renewable energy, representing between 70% and 80% of scientific production [15]. In general, countries are in the process of achieving, and in some cases exceeding, many of the objectives set in their commitments under the Paris Agreement. This is enough to reduce the expected increase in global CO2 emissions related to energy, but it is not enough to limit the global warming to less than 2 ◦C. The increase in energy consumption, especially in HVAC systems (heating, ventilation and air conditioning), and CO2 emissions have contributed to more energy policies focusing on the development of energy saving and efficiency strategies as an absolute priority [16]. An example of this is the European Energy Performance of Buildings Directive.

The major advanced economies: USA, the European Union, and Japan appear to be clearly on track to achieve their climate commitments, although it will be vital that these countries introduce other additional improvements in energy efficiency. A clear example of it would be the cement production industry, one of the most polluting. However, several mitigation measures are gaining importance in the recent years, and this industry is increasingly turning towards cleaner production. This is because new alternatives are used as raw materials and fuels. In existing industries, mathematical statistical models, simulation processes, optimization of procedures are everyday methods [17]. Natural gas and oil are going to continue to be the essence of the global energy system for years to come, as governments have indicated in terms of climate. Faced with this situation, the fossil fuel industry should consider the risks of a more direct transition. The interdependence between energy and water will become more pronounced in the upcoming years. In this sense, the management of the relationship between energy and water is crucial for the successful fulfilment of a series of development and climate objectives. It is evident that, to promote energy saving, it should be tried to create efficient energy systems and banish old ones. In fact, there is great potential for energy savings in the renovation of the housing stock in many cities, as is the case in Denmark, where 75% of the buildings were built before the eighties, when the first important demands were introduced for the energy performance of the building [18].

In this way, without altering the development of the technique, significant energy savings will be achieved. Nor should we ignore the human factor; the education and global awareness of citizens is important to achieve energy savings.

Recent studies show how construction technologies have been developed to improve comfort and energy savings in both conventional buildings using concrete or heavy structures, and also in timber buildings. An example of these technologies is the thermal insulation added to the facades and the improvement of the performance of the glazing or a district heating designed in a residential area. Nevertheless, although these technological advances can be implemented in concrete and wood buildings, the studies carried out highlight how timber construction is associated with thermal comfort and is perceived as innovative [19].

This is paradoxical because, after having expanded the use of gravel and concrete, the wood construction was increasingly neglected for a period of time. However, in recent years the Kyoto Protocol has pushed to limit CO2 emissions in all productions. The use of sustainable materials is definitely being stimulated worldwide, as they can store carbon dioxide. The advantages of using sustainable materials have been confirmed by various scientific studies [20].

Buildings using sustainable materials are called green buildings. They reduce and optimize their energy and water consumption as much as possible and integrate into their environment, whether natural or urban, causing the least possible environmental impact. Smaller buildings tend to achieve these objectives satisfactorily. However, in larger green buildings it is much more difficult to meet the requirements of comfort, lighting, noise, use of space, or more hours of operation. It is feared that these innovative buildings will run the risk of repeating past mistakes, especially if they are too difficult to manage [21].

In this process of constructing eco-friendly buildings, it is fundamentally important to know the user's perception of indoor environmental quality and satisfaction in green buildings. Numerous studies show that employees perceive green spaces better than conventional buildings, but those surveyed are not able to differentiate between both types of buildings. However, there is an unequivocal correlation between ecological attitudes among employees of green buildings, showing them to be more satisfied and more prone to sustainable attitudes [22].

In order to achieve appropriate low-carbon rehabilitation interventions, users play a very important role. A sound learning process must be ensured for occupants, designers and building managers. The behavioural habits of citizens are themselves an essential component in reducing consumption patterns [23].

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

The aim of this work is to analyse all the world scientific publications of renewable energy indexed in the Scopus database in order to study the most important areas of investigation in energy saving, as well as to highlight the existing worldwide relationships between researchers and institutions. It is decided to select Elsevier's Scopus for this study as it is the largest database in the world. Thanks to the API interface, it is very easy to automate the search for literature related to the topic by selecting it by authors or by institutions. The two main scientific databases, Scopus and Web of Science pose the question of the comparability and consistency of the results of the different data sets [24]. With respect to journals reported in the two main scientific databases Scopus and Web of Science, a comparative study reveals that the volume of journals on the Web of Science is inferior to the number of journals in Scopus [25], and the relationships obtained in the results with both databanks for the items and the citation numbers by country, and for their rankings, are very robust (R2 ≈ 0.99) [26]. This methodology was successfully used in several scientific fields, so only Scopus data has been used in this research [27].

#### *2.1. Automated Data Collection Applying Scripts*

The flowchart shown in Figure 1 explains how automatic data extraction would be carried out from Scopus. The process that is carried out by Research Network Bot (ResNetBot) [28] could be separated into the following stages: (1) Obtaining information on documents containing 'energy savings' as keywords or included in the title. We must notice that analyse the keywords in scientific articles is one of the most common research topic in different fields of engineering science. This information is represented graphically, where the size of the nodes is proportional to the h index and the lines between two nodes indicate that they are cited. (2) Obtaining information from authors of the works mentioned in step (1), which includes for each author the Scopus author identification number, affiliations, publications and dates, number of citations and h-index. (3) Processed to obtain the collaborative network of authors who have published articles that contain 'energy savings' as keywords or are included in the title. The information stored refers to the number of collaborations between co-authors and their affiliations, country and city. With all this information a graph will be created where the size of the nodes will indicate the scientific production according to the country and will be proportional to the H index and the distance between them will represent the collaboration with one another.

**Figure 1.** Flow diagram showing how the software obtains information from the Scopus database (ResNetBot).

#### *2.2. Data Evaluation and Polishing of Texts*

ResNetBot obtained data that was classified by domains and stored as a compendium of text files in JSON format., http://www.json.org. The Scopus API gives the option of requesting information through all type of details, so we programmed the bot to request all the information available in the files and choose the most convenient data for our work. As a result of the vast volume of data, some irregularities were observed during the work of checking the information. This, however, is a frequent and common problem in these large databases [29]. Phrases with a very similar meaning may sometimes appear written with slight changes. For this reason, it is very important to carry out a debugging process that allows us to obtain the appropriate information. The same expression may not be exactly the same written depending on whether we use capital letters in the first word, the second word, or both. We see an example of this, "Energy saving", "energy Saving", and "Energy Saving".

Therefore, some of the improvement criteria are applied from OpenRefine (http://openrefine.org), which contains some formulae such as clash of keys and neighbourhood algorithms to combine all expressions that correspond to the same idea but are not written in the same way. In order to classify data information and identify singular values, data sheets are finally used; this program has

already been successfully carried out in former studies [30] and could be performed on key-words and author names.

#### Diagnostics and Data Viewing

The data collected by ResNetBot and further processed with OpenRefine was archived in a separate databank. Once this process has been completed, the latter information will be analysed using visualisation tools and statistics based on diagrams.

A graph is a set of vertices(elements) joined by edges. Graphically, the vertices are points and the edges are lines that join them, representing a relationship between two vertices. The visualization of graphs is a great help to determine the interrelation between the different parts. On the other hand, this representation makes it possible to incorporate specific values of the elements and their links. Over the last few years, powerful graphic representation programs have been created that enable the exhaustive analysis of the particularities of the graphics through the use of various functions, such as the reorganization of the vertices and the design of the vertices represented with different colours and sizes depending on their properties. One of the best known when using open source instruments is Gephi [31], which encompasses several statistical tools.

#### **3. Results**

#### *3.1. Energy Saving Publication Relationship*

In Figure 2, the documents obtained in the search described above have been represented all together. The graph represents the relationship between the documents, with each node representing a file. Of the total of 20,095 papers, in the outer halo formed by 12,402 individual papers are represented those that do not cite any other Energy Saving paper, this is 61.7%. Then there is another intermediate halo of 1919 papers that are linked to at least one other paper, but do not connect with the central core, which represent 9.6%. Finally, there is the central nucleus formed by 5774 nodes that are linked together, and that account for 28.7% of the total, is what could be defined as works that scientifically deepen in that field, because they are based on the previous ones to make some type of contribution.

In our research, there are a total of 70,951 keywords. As there is a large number of papers that are not part of the central core, we proceed to a refining of the 5774 core papers. 23,083 unique keywords were obtained, which after the refining process, 10,119 different keywords remain. After this process, the analysis of the documents that are related to each other is done, and they are represented according to the communities detected with the modularity algorithm of Gephi. Once this process is finished, Figure 3 is obtained. Of the 29 communities detected, the 15 main ones, with a percentage higher than 2%, have been represented in colour to distinguish them.

**Figure 2.** Display of the whole energy saving publication relationship indexed by Scopus.

This process has been necessary to reduce the excess of information that would impede the global vision we pursue.

Figure 4 shows the main 15 communities which have been named to facilitate the analysis. This table shows the main keywords of each community detected.

**Figure 3.** Display of the whole energy saving publication relationship indexed by Scopus.

Fifteen communities are detected within the large network of publications related to the term energy saving. In Figure 5, it can be seen that the area with the highest number of publications is lightning and thermal comfort, followed by emission reduction and building insulation. The high energy consumption that occurs in homes and businesses, together with the growing concern about the impact that energy consumption has on the environment, could show how important the search for facilities and systems that enable energy savings in these field is becoming. On the contrary, the railway and public buildings communities are in the last place of the ranking, with a percentage lower than 4%.

#### *Energies* **2019**, *12*, 2030


**Figure 4.** Communities detected for Energy Saving.

**Figure 5.** Communities are detected within the large network of publications related to the term energy saving.

#### *3.2. Authors, A*ffi*liations and Countries in Energy Saving*

Regarding the authors and their relationship among them, it has been seen that in total there are 43,860 researchers. In Figure 6, the authors have been represented by colours according to the country of origin. This figure shows that the researchers that contribute most to the development of scientific production, come from academic and scientific institutions from China, which has the highest number of authors with more than half of the total number, followed by Italy (9.5%), the USA (8.9%), and Japan with 6.8%. Afterwards, with less percentage are Malaysia (4.77%) and Spain (2.2%), and finally a block of countries with a percentage below 2%. Null means that the Scopus register does not contain information on the country of the author.

**Figure 6.** Collaborations between scholars from various countries. The size of the nodes is proportional to the H index.

From Figure 6 it can be drawn as conclusions that among Chinese researchers there is a tendency to collaborate among themselves, but also at the international level with American authors, Malaysians, and to a lesser extent, Japanese authors.

Whereas at European level, the relationship between its authors is with those authors belonging mainly to the European environment. These appear isolated, as can be observed in the lower left part of the graph, where it is extracted that the Italian researchers, the most productive, are mainly connected with the Spanish or with Greek authors among others, but also Japanese and Canadian.

Several layers are observed, being the most external (as explained above) that of the authors who have written about energy saving, but who do not collaborate with anyone in energy saving. Here are 2763 authors, which means a large minority. Then there are intermediate layers, where groups of collaborations are going from less to greater intensity. In this intermediate layer there are 41,097 authors, and it is observed that it is the most numerous. As a result of this international collaboration, many papers have been published.

The method is focused not only on the analysis of the scientific collaborations between institutions of each country, but also to determine which of these institutions are more involved in energy saving.

Researchers associated with Chinese institutions is shown in Figure 7, highlighting the number of people involved in energy savings in the following order: Tsinghua University (2.27%), North China Electric Power University (1.98%), and Tianjin University (1.91%).


**Figure 7.** Connection of authors belonging to Chinese institutions.

In Figure 8 and Table 1 we can see the evolution of the publications developed by the most productive institutions worldwide in the last 10 years (2009–2018). On the bottom of the figure, the production of Ministry of Education of China can be seen, which is the most important institution registered. It does not surprise that, as a fact, it would englobe a lot of researching in the most populated country in the world. This institution represents the 1.64% of the total research (718 publications). Above it, the Tsinghua University represents the 1.59% (699 publications). This institution is one of the biggest universities in Asia with more than 40,000 students, 32,000 of which belong to masters or doctorate degree. After them, we have more institutions and as we can guess, all of them belong to China. It was a highlight surprise to notice all the research developed in this country and the high difference between China and other countries, at least in this search.

**Figure 8.** Ten most productive institutions during the period 2009–2018.


**Table 1.** Number of publications of the ten most productive institutions and universities between 2009–2018.

*3.3. Typology of Literature, Languages of Publications and Distribution of Output in Subject Categories*

Table 2 shows the different types of publications during period 2009–2018. We can see how more than half of the total of publications of the search are articles (50.7%). Following them, the conference papers, with 43.1% of the total. After that, the quantity decreases to different types of documents as Reviews (2.3%), Book Chapter (1.6%), or Conference Revision (1.05%).

**Table 2.** Distribution of the searching documents according to the type of publications between 2009–2018.


The vast majority of articles are published in international journals. Because of that reason, more than 87% of the work have been published in English (31,992), approximately 10% have been published in Chinese (3642) and only 3% includes other languages as Japanese, German, Russian, Polish, or Korean. In Table 3, shows the publications organized by language. Figure 9 shows the classification of publication results by subject, according Scopus.


**Table 3.** Most used languages in documents.

**Figure 9.** Classification of publication results by subject, according Scopus.

#### *3.4. Literature Distribution by Area and Entity*

Figure 10 shows the distribution of publications by country. On the map, red is associated with a high number of publications. Green indicates that there are less publications or grey means that there are no publications or a very limited number of them. It is easily observed how China is the country with more publications (18,309), followed by United States (3716), Japan, and Italy. Also highlight Germany, United Kingdom or Taiwan. This information leads us to the idea of how important the energy saving in the industrialised countries has become, as the richer states have also become involved in this issue by using harsh policies to achieve the proposed objectives. Figure 11 shows the 15 countries that have contributed the most in the last decade.

**Figure 10.** World Map with the 25 countries with a higher number of publications.

It should be taken into account the number of inhabitants of each country in order to properly analyse the results. It is not a big surprise to find out that China is the country with more publications as the population in China is the highest in the world (1,386,000,000 in 2017). This is followed by United States, also with a high population: 327,200,000 at 2018 and Japan with 126,800,000 in 2017. However, it is interesting to find countries such as Italy (60,590,000 in 2017) or Germany (82,790,000) with a lower population but still with a significant number of publications.

In Figure 12, the five keywords that appeared the most in our search are shown. 4/5 of these are related to energy, being the usual investigation and research of this field to create and find new ways of energy optimization and its use in new technologies. We can see how the term "Energy Conservation" started to be used in the beginning of the century and it has increased remarkably from 2013. In the second place, there is the term "Energy utilization" followed by "Energy Efficiency". After them "Energy saving" was very often used from 2008–2010, but it decreased as the researchers tend to focus their work in the conservation of the energy more than the saving of it.

**Figure 11.** 15 countries that have contributed the most in the last decade.

**Figure 12.** 5 most important Keywords and how often they were used during the last decade.

In Figure 13 and Table 4, the ranking of the six most important journals are shown. We have analysed not only the Scimago Journal Rank by Elsevier, but also the CiteScore. The metric of the latter calculates the citations of all papers submitted in the past few years preceding a title. It also includes Source Normalized Impact per Paper, SCImago Journal Rank, quotation and documents counting, and quotation rates. In the figure we have included also the H-Index as it has become very popular in the last decade and the reader can figure it out by the size of the circle according to each of the journals. We can see how "Renew. Sust. Energ. Rev." is the journal with higher not only SJR (3.04), but also CiteScore (10.54).

**Figure 13.** Classification of the first 6 journals. CiteScore is the Scopus ranking and SJR (Scimago Journal Rank) by Elsevier. The size of the circles indicates the H-index.


**Table 4.** Most important journals through the research.

The 21st century is the century of sustainability. This is why the building sector is investing a lot of resources in increasing energy efficiency, concerning both the materials and technologies used to make buildings and building techniques more energy efficient. In other words, it is not only a question of constructing more energy-efficient buildings, such as passivehause buildings and internship rehabilitation, but also of how to construct them, applying more efficient construction techniques, energy savings that must also be extended to the manufacture of construction materials and their transport to the construction site.

In Figure 14 we can see the most important authors related with the search through 2009 to 2018. The most important author is Kansha, who against all odds, is not from China, he is from Japan and has published a total of 57 papers related to the search, representing the 0.16% of the total of publications. He has a total of 1103 citations through all his work and a h-Index of 21. After him, Yang has a total of 47 works (0.13%) related to energy saving. All his work has a total of 8829 citations, getting a 45 H-Index. Other authors such as Li go after with 45 works (0.12%) and the Japanese Tsutsumi with 44 (0.12%).

**Figure 14.** Authors who have published related with the search through the last 10 years.

#### **4. Discussion**

The scientific publications in relation to energy saving have had an exponential growth since the year 2000, from that date until now more than 90% of scientific production is condensed. This indicates that the last 20 years are those that have generated the interest of the States through their researchers and programs to support this research. For example, the research related to energy saving in communication and computation [32]. The works published in the topic of energy saving, are mainly in the form of scientific articles (50.7%), but there is an important number in the form of conference papers (43.1%), which indicates that today it is a novel field of research, in which this topic is addressed in many scientific meetings. It is observed that the books are a percentage below 1%, which means that it is not yet a very well established issue.

Regarding the language of publication, as expected, is dominated by English in more than 87.1%. However, when compared to other scientific fields, it is usual to be above 90%, which suggests that journals in particular countries address this issue in a significant way. As example the Dianli Xitong Zidonghua/Automation of Electric Power Systems, or Zhongguo Dianji Gongcheng Xuebao/Proceedings of the Chinese Society of Electrical Engineering. Such is the case of China, which contributes almost 10% of all publications in its own language. To a lower extent, German (e.g., the journal eb—Elektrische Bahnen) and Japanese (e.g., the journal IEEJ Transactions on Power and Energy from The Institute of Electrical Engineers of Japan), both around 2%, stand out. In any case, it is clear that China, as one of the world largest consumers of energy, has a great concern for energy saving [33], as was to be expected. This is also evident if it is observed that the main institutions that publish on this subject are also from China, Table 2. If the analysis focuses on the scientific fields where these works are assigned, Figure 9, it is remarkable that the first is not energy, but that the first place is occupied by engineering with 33.1%, showing that developments in this field are the largest part of the contributions, especially those related to energy conservation [34]. On the other hand, it is surprising that Computer Science plays such an important role in this field; this is justified by optimisation-related works [35]. Finally, it should be noted that one of the largest group is "others" with 9.2%. This should indicate the great diversity of scientific fields from which this subject is approached, being important the contribution to the energy saving that is studied from the greater number of points of view [36], which is what will allow to contribute to the global energy saving including geographical areas with serious energy shortages [37].

From the point of view of the networking of the works among themselves, this work shows a vast number of works published on the subject of energy saving, but the most important fact is that about 62% do not cite other work related to this subject. This implies great individualism, and scarce scientific connection in this field, which should make researchers reconsider, because it is

important that research is based on previous studies to progress. However, it may also imply that they are researchers from other fields of knowledge who are becoming aware of or are including these items in their work, which is very encouraging. Within the works that are related to each other and are the core of scientific research in energy saving through the identification of communities have been detected the top 15. In the central positions, we find the scientific fields that might be said are most popular, see Figure 4: Lighting, thermal comfort, building insulation, cooling techniques, HVAC techniques, household-behaviour, or public buildings. Lighting has undergone a revolution since the emergence of solid-state light sources. The efficiency of solid-state technology provides, in a multitude of applications, energy savings, and environmental benefits [38]. Above all energy saving based on lighting has been studied for residential issues [39], public buildings, schools [40], or shopping malls [41]. Another one of the great challenges of energy saving is thermal comfort. Thermal comfort implies a good indoor climate that is important for the success of a building [42], not only by making its occupants feel comfortable, but also by deciding their energy consumption and thus influencing its sustainability [43]. On the other hand, persons have a natural trend to adapting to the conditions changing in their surroundings. This is especially the case with the scientific communities found around the central core of energy saving: building insulation, cooling techniques or HVAC techniques. Therefore, saving energy from households is one of the priorities of most countries, such us China, Italy, USA, or Japan. Those that occupy external positions are considered less integrated with the others, still have great possibilities of integration with the core and are generally the smallest. This is the way to find it: Distillation, communication protocols, electricity market, communication equipment, or railway. Related to these scientific communities it should be mentioned that distillation is still the most popular separation technology of chemical plants, and it uses around half of the operational costs, which shows the major concern of the researchers in applying here the energy saving. The energy savings in the communication protocols has been developed mainly using mobile telephony, where the terminals have a limited battery capacity. The third scientific community found is that related to the electricity market, being this one of variable tariffs, makes energy savings in certain time slots is of special interest. Regarding communication equipment, given the rising digital level of the developed countries, there is a growing concern about the data centres, especially in two respects, cooling and IT load. The amount of world energy consumption has gone from 0.5% in 2000 to 1% in 2005. With respect to the railway, the electric railway system is considered the most environmentally friendly means of locomotion. This is due to two important points: Firstly, the low resistance to running and secondly, the regeneration of energy during braking.

#### **5. Conclusions**

This study presents a wide range of data on the international contribution to scientific knowledge in the field of energy savings during the period from 1939 to 2018. The newness of this work is the application of community detection algorithms which identify scientific communities in the world that develop their work in fields related to energy saving. They also detect the collaboration preferences of their authors at an international level or even at an institutional level as in the case of China. It should be noted that many of the published papers have been the result of this international collaboration.

This methodology makes it possible to understand two complex systems: the network of publications based on quotations and the network of scientific collaborations based on the nationality of the co-authors.

Research Network Bot was used to search for publications indexed by Scopus containing 'energy savings', in the title or keywords. The information obtained was used for the analysis and graphical representation of these networks of scientific collaboration of researchers working in this field.

49,539 documents have been found, from which 55.5% are articles and 38.5% are conference papers. In this period there has been an exponential increase in publications, highlighting the categories of Engineering, Energy, and Computer science.

In light of the data obtained, Chinese institutions have provided an important contribution to the field. It is noted that, logically, the compound terms related with energy are the most used keywords, being energy conservation, energy utilization, energy efficiency, energy saving, and fifthly, optimization, proved as the most popular terms. Different aspects of the publications are analysed, such as language, type of publication, scope, top contributors' countries, institutions, journals, authors, as well as the frequency of appearance of keywords. Geographically, the giants China and United States, followed by Japan and Italy are highlighted as the top contributors' countries. The most active categories within the field of Energy are Lighting with 9.44%, Thermal Comfort (9.18%), Emission reduction and Buildings Insulation, both of then representing more than 8.5%. These 4 fields represent areas where energy saving plays a very important role as buildings, responsible for more than 40% of the energy consumption. In order to reduce this, there should be a release of a kind of regulation or policy to control systems for the installation of lighting.

Therefore, the communities that are in the lead, refer to the state of comfort, but also to the latent concern for the waste emitted into the atmosphere and how these can directly affect our planet in the future. Policies actually adopted in developed countries are insufficient to ensure a significant reduction of the energy consumption. It forecasts that the current global demand for energy in buildings will double by 2050. Thus, it is necessary to implement, at the global level, effective saving policies [44]. The challenge is always to reduce the energy consumption without forgetting the comfort of the user and the requirements of the building. The results of the work show that during the last 10 years have substantially increased the number of publications in this field, which indicates that research in terms of energy saving is assumed as an important and relevant issue in the international scientific scene.

During the last decade China, the most productive country in "energy saving" has become the main leader in the growth of energy efficiency. Its GDP grew a total of 7.8% per year and the energy intensity fell a total of 18.2%. This outstanding result is related with how energy efficiency policies have become the center of the political landscape through these year [45]. Energy efficiency was first included in the five-year plans, in the mid-2000s. It was a first sign of China's acceleration towards greener policies. In addition, energy efficiency has played a role in number of key Chinese economic policies: The market liberalisation and privatisations of the 1990s, the quest for energy security, economic development, or the fight against climate change.

This shows that there is a correlation between the number of publications in this country and the application of energy policies, but we do not have data to extend this assertion to the rest of the countries of the world. The study of the implementation of energy policies in each country, the energy savings achieved and the relationship with its scientific production on the subject, both quantitatively and qualitatively, would be the subject of future work.

**Author Contributions:** J.L.d.l.C.-F. and C.d.l.C.-L. performed the literature review and article writing. A.A. and A.-J.P.-M. analysed the data using community detection methods and article writing. F.M.-A. and F.G.M.: Research idea, article writing, and formatting; they jointly contributed to the structure and aims of the manuscript, paper drafting, editing, and review. All authors have read and approved the final manuscript.

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

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

#### **References**


© 2019 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*

## **Factors That Contribute to Changes in Local or Municipal GHG Emissions: A Framework Derived from a Systematic Literature Review**

#### **Isabel Azevedo \* and Vítor Leal**

Institute of Science and Innovation in Mechanical and Industrial Engineering (INEGI), Faculty of Engineering of the University of Porto (FEUP), Rua Dr. Roberto Frias, 400, 4200-465 Porto, Portugal; vleal@fe.up.pt **\*** Correspondence: iazevedo@inegi.up.pt; Tel.: +351-229578710

Received: 15 April 2020; Accepted: 16 June 2020; Published: 19 June 2020

**Abstract:** The changes observed in a municipal or local energy system over a period of time is due a number of concurrent dynamics such as the local social and economic trends, higher-level policies (e.g., national), and the local policies. Thus, a thorough identification and characterization of the many factors of change is key for an adequate assessment of the effectiveness of policies adopted at the local level, as well as for future planning. Such thorough identification and characterization of the factors that are determinant for the evolution of the local energy system is the content of this paper. The identification is performed through a systematic literature review, followed by a synthesis process. The factors identified are grouped into the categories of local context, local socio-economic and cultural evolution, higher-level governance framework, and local climate change mitigation actions. They are represented through a set of measurable variables. The factors and respective variables can be used to improve the disaggregation of changes in the local energy system into individual causes, leading to a better assessment of the evolution over time.

**Keywords:** local climate policy; municipal authorities; municipal energy planning; local GHG emissions; multilevel governance

#### **1. Introduction**

Local actions are crucial for the achievement of global climate change mitigation goals. Both the physical characteristics of the local energy systems and the regulatory competences of local authorities prompt local level as an appropriate level for action [1]. Moreover, local actions are becoming more common, due to the increasingly active role of local authorities and to the recognition of their importance by policy makers at higher governance levels. Thus, there is an increasing need to assess their contribution towards climate change mitigation, and to the achievement of global targets.

However, given the universal scope of climate policy and the link between the local energy system and local context, it is currently very difficult to properly assess the actual contribution of local policy actions. Indeed, the local energy system is influenced by numerous factors, of which only some are related with actions promoted by local authorities. The changes in the economic activity or on population are two examples. The assessment difficulties are amplified by the fact that local policies and actions are implemented as a part of a complex multi-level governance framework.

Previous studies have already highlighted the importance of considering external factors when evaluating local actions. Within the studies reviewed in [2], several factors were identified as key in the link between observed changes and actual effects of local climate change mitigation actions. For instance, the evolution of local socio-economic features (structure of local economy, demographic and economic conditions of local population, etc.), prompted by economic, social, and technological changes, is one of the aspects that are commonly considered relevant for the observed variation in local GHG emissions.

However, within existing studies, a comprehensive analysis of all the relevant factors that influence the local energy system and the corresponding GHG emissions was not found.

Moreover, the analysis of empirical studies focusing on assessing the effects of local actions has emphasized this need to thoroughly identify all the factors that contribute to changes in local GHG emissions [3–5]. For instance, in [3], Millard-Ball performed a quantitative analysis on the effects of city climate policy, without considering any external factors; although some scattered evidence of climate plans effects was found, the findings are not sufficiently robust (which could be partly explained by the non-consideration of other relevant factors). The evidence for impact of local climate plans in the different sectors is not demonstrated with the empirical analysis, which does not provide a coherent result across the different models that were tested. Similarly, in [5], the conclusions drawn were also insufficient to prove the causal link between local actions and GHG emissions. Within the results from the linear regression analysis with panel data, the coefficient associated with the effect of local actions is not statistically significant for any of the countries that were included in the study. However, in the same analysis, the relevance of municipality specificities as well as of higher-level context (i.e., country) for the assessment of local GHG emissions and their evolution over time was confirmed (resulting in coefficients that are statistically significant). Both these studies show the need for a comprehensive consideration of the different factors that impact local GHG emissions to ensure a robust and meaningful assessment of the specific impact of local climate change mitigation actions.

Thus, the identification of the factors of change of local GHG emissions is a crucial step towards a more complete understanding of the observed changes within the local energy system and assessment of the actual effects of local climate change mitigation actions. Indeed, to develop a model that is able to adequately represent the evolution of the local energy system, identifying the links between the changes in local GHG emissions and the different explanatory factors, it is necessary to first identify all the factors that may have an impact on local GHG emissions (factors of change). Herein, factors of change refer to all features and events that lead to changes in the local energy system and respective GHG emissions. These include factors related with local climate change mitigation actions as well as external factors (i.e., factors that are not directly influenced by actions promoted by local authorities). For instance, local demographic evolution—size and structure of local population—may lead to significant changes in the local energy system and, consequently, on the respective GHG emissions.

This paper attempts to perform a systematic identification and characterization of the factors of change of local GHG emissions. This comprises the identification of all factors that may potentially influence local emissions, even if their specific effects in the local energy system are not quantifiable or even relevant (depending on the local system being studied). In order for this to be applicable in the future, it is necessary not only to identify the relevant factors of change, but also to characterize them with indicators. Therefore, that is also covered by this paper. The work is presented in five sections. Section 2 describes the methodological approach that was used to identify the factors of change of local GHG emissions and the review process, including the documents used as reference to this analytical framework, as well as the presentation of the final list of factors identified. Section 3 is dedicated to the presentation of the results, including the description of the identified factors of change and the listing of the respective quantifiable indicators. Section 4 presents a brief discussion of the obtained results, focusing on their potential contribution for improved planning and evaluation of local climate change mitigation actions. The paper ends with some final remarks regarding the results obtained (Section 5).

#### **2. Methodology**

Considering the extensive research that has already been done on this topic, a wide literature review was considered an adequate methodological approach to identify and characterize the factors of change of local GHG emissions. More specifically, the method used consisted in a systematic literature review, using content analysis—which allowed reviewing, assessing, and aggregating a body of literature utilizing pre-specified standardized techniques to all the documents reviewed [6]. The definition of the specifications for data collection and analysis beforehand and respective use as a guide for performing the process aimed to avoid bias in the review.

A systematic literature review process consists of these four sequential steps: (i) definition of the objective; (ii) definition of the specifications to be followed; (iii) the retrieval and assessment of the information; and (iv) synthesis of the results obtained. In what concerns the specifications defined a priori as guidelines for the review, these consist in three types of rules/criteria: (1) the eligibility criteria to retrieve the literature; (2) the rules to retrieve information from literature; and (3) the rules for the analysis and screening of the obtained information. Table 1 presents the details of these three types of criteria.


**Table 1.** Systematic literature review specifications.

The selection of relevant research covered the following keywords: "Energy system characterization"; "Energy system modelling"; "Energy system planning"; "Energy system scenarios"; "Energy policy evaluation"; and "Local climate change mitigation". The keyword selection was based on the different contexts within which this discussion has been raised, including planning, policy evaluation, and energy management.

To guarantee the reliability of the results, the documents considered in the review were restricted to peer-reviewed articles published in scientific journals or conference journals, reports or working papers published by internationally reputed entities (as e.g., International Energy Agency and European Union), and PhD theses. The reviewed documents were published between 1999 and 2017. A list of the documents with a description of their main focus of research is presented in Table 2.

**Table 2.** Characterization of the documents used for the identification of the factors of change.



**Table 2.** *Cont.*


**Table 2.** *Cont.*

Legend: JP, peer-reviewed journal article; WP, working paper; CP, conference paper; Bk, Book; PhD, PhD thesis; PR, policy report; TR, technical report.

The information collected from the documents comprised three types of relevant information: (1) factors of change of GHG emissions; (2) variables that may be helpful to characterize the factors of change identified; and (3) information on the relation between the different factors of change and local GHG emissions.

In regard to the identification of the factors of change, Figure 1 presents a schematic representation of the whole process. The analysis of these documents resulted in a list composed of 511 variables with an impact on GHG emissions. The full list, as well as the reference to the documents where each variable has been referred, is presented in Table A1 (Appendix A).

Then, to reduce the list to the essential, both the eligibility criteria (Box 2 in Figure 1) and the screening filters (Box 3 in Figure 1) were applied successively. The following paragraphs describe this process in detail.

The eligibility criteria were defined to guarantee the applicability of the framework to the assessment of local climate change mitigation actions. Firstly, it was assumed that the respective change must be associated with energy-related GHG emissions, excluding, e.g., the emissions from crops and livestock. Then, the associated effect on the energy system should also be on Scope 1 and 2 emissions—GHG emissions that can be attributed to each municipality (i.e., that local inhabitants can be made accountable for). For instance, with this condition, changes in aviation and international commerce emissions are excluded, since typically an airport serves an area much larger than the municipality. Finally, it was also defined that factors should either change over time and/or influence the effect of the other factors identified. As the goal is to decouple the observed changes in a specific local energy system into effects of the different factors of change, it makes sense that factors that are constant over time (even if determinant for local energy use) must be disregarded. For example, topography of the local territory can have significant implications for local energy use, but it cannot be considered a factor of change as it cannot lead to changes of a specific local energy system over time.

**Figure 1.** Schematic representation of the systematic literature review process.

The implementation of this step corresponded to the exclusion of variables that do not comply with at least one of the three eligibility criteria. The first eligibility criteria led to the exclusion of variables such as the "existence of strategies for landfill methane capture" and "changes in cattle production or fertilization strategies". Moreover, the constraint that implies that each factor must influence local GHG emissions led to the removal of variables such as "level of international trade" and "trade allowances between regions". Lastly, the criteria that exclude variables that do not change over time nor have an impact on policy outcomes implied the exclusion of factors "latitude", "topography", and others. Overall, with the application of the eligibility criteria, 73 variables were excluded from the list. The full list of the non-eligible variables is presented in the Appendix A (Table A2).

Furthermore, a screening process was defined and implemented to reduce the list of identified factors to the essential. This screening process is composed by three main steps that correspond to the application of three distinct filters: (1) explicit repetition; (2) implicit repetition; and (3) redundancy. The screening of the list of variables led to the exclusion of 331 variables. The full list of the excluded variables and the identification of the filters that led to their exclusion are presented in the Appendix A (Tables A3–A5).

The filter of explicit repetition refers to the removal of duplicates, i.e., factors that were identified by more than one document. For example, the number of inhabitants (or overall local population) and GDP were two of the variables that appeared more than once in the list.

Implicit repetition refers to factors that, even if not literally in duplicate, if withdrawn from the list do not imply an information loss. This includes factors that are a combination of other two or more factors that are also listed. For instance, average household area per capita can be considered as implicit repetition if both overall population and built area for residential purposes are already listed. This led to the exclusion of variables such as average household area per person (as total household area and overall population were also listed) and sectorial GDP (as GDP per subsector was also included).

Finally, the redundancy filter corresponds to the screening for dependent factors, whose relation has already been identified by literature. With this filter, for instance, average income per capita was excluded, considering that changes in private consumption (also linked to average income) would imply the exact same impact on local GHG emissions.

Once the variables with an impact on local GHG emissions have been identified and screened, a final list that is consistent with the goal of this work was obtained. The application of the eligibility criteria and the screening filters led to a final list of 107 variables that were grouped into 40 factors of change. This was done through an aggregation process, distinguishing policy-based factors from the ones that occur naturally, and factors related to actions implemented at local level from factors associated with governance actions at higher-levels. This grouping process is detailed in Table 3.


**Table 3.** Final list of variables and respective grouping into factors of change.


**Table 3.** *Cont.*


**Table 3.** *Cont.*

Finally, to complete the theoretical framework, each individual factor was characterized with (preferably quantifiable) indicators. For most of the factors, the documents where they were extracted from also mentioned possible indicators. When more than one indicator was available, the choice fell back on the indicators whose data are more commonly available. When the documents did not mention possible indicators, these were searched for in national statistical agencies and policy reports from national or international governments. Within this step, there was an effort to choose indicators that use data that are commonly available for the municipality level, or equivalent scale, to avoid the need for downscaling.

#### **3. Factors of Change of Local or Municipal GHG Emissions and Respective Characterization**

This section presents the results obtained by the systematic literature review process described in the previous section. This includes discussion detailed description and an analysis of the final list of factors of change that were identified as well as the listing of the characterizing indicators for each factor of change.

To better understand the relation between the different factors of change and prompt a more fulfilling discussion, the 40 factors of change identified in the review were grouped into four distinct sets. Three of these groups correspond to three sources of change: local socio-economic and cultural changes; higher-level governance framework; and local climate change mitigation actions. The fourth group (herein referred to as local context) comprises variables that, despite not leading directly to changes of local GHG emissions, have influence on how the different sources of change impact local GHG emissions (e.g., climate awareness). A schematic representation of the relation between these sources of change and the observed changes in local GHG emissions, considering the local context, is presented in Figure 2.

**Figure 2.** Framework of the sources of change and the observed changes in local GHG emissions.

#### *3.1. Local Context*

Local context refers to the local characteristics that, despite not leading directly to changes on the local energy system, influence the outcomes of both local and higher-level policies. On the one hand, this comprises the local inhabitants' characteristics, as education level, cultural diversity, and awareness towards the problematic of climate change. On the other hand, local context also includes characteristics referring to the local authority as an entity, including human and financial resources as well as legal competences. In general, the characteristics here referred as local context do not change significantly in the medium-term.

A complete list of the factors of change considered as local context is presented in Table 4, along with the respective characterizing indicators.




**Table 4.** *Cont.*

A summarized description of each of the factors of change associated with local context is presented below.

#### 3.1.1. Local Inhabitants' Characteristics

Education level refers to the average level of education of the local population, e.g., whether the majority of local citizens have completed high school, or even a bachelor degree. Here, it was assumed that the distribution of local population per education level could be a good indicator.

Cultural diversity refers to the heterogeneity of local population in terms of cultural background. This is mostly linked with the scale of immigration to the city, currently but also in the last decades. Culture is associated with daily habits and, thus, with specific patterns of energy use and reaction towards certain policy actions/incentives. One of the possibilities to assess this diversity could be through the number of residence permits by reason (or by country of origin).

Willingness to act refers to the level of motivation of the local citizens in having an active role towards climate change mitigation, which is highly associated with their level of awareness to the problem and their civic involvement in the community. For instance, the level of awareness of local population on environmental and climate change issues can have a huge influence on the outcome of any related policy independently on the level of implementation. Even if waste associated emissions are not considered within this work, it is assumed that the share of waste separation could be a good indicator of the local population awareness on climate change issues. Moreover, the degree of civic involvement may be assessed by the number and scale of environmental nonprofits that were created at the local level.

#### 3.1.2. Local Authorities' Characteristics

Legal competences refers to the regulatory competences of local authorities regarding climate change mitigation, i.e., their scope of action. These may vary significantly from country to country and are determinant for the type of actions that are promoted at the local level. The level of autonomy of local authorities can be seen as a good indicator of their scope of action.

Qualified personnel corresponds to the technical competences of the municipal staff for developing and implementing actions targeting climate change mitigation. This can be assessed through the existence of a dedicated agency or municipal department to climate change mitigation, and from past experiences with these issues.

Economic situation translates the economic capacity of the local authority in financing and promoting actions towards climate change mitigation. This can be assessed by the relation between income and expenses of the local authority, which provides a characterization of their yearly profit or debt. Moreover, the ability to obtain funding, through private partnerships or public tenders, is also a relevant indicator.

#### *3.2. Local Socio-Economic and Cultural Changes*

Local socio-economic and cultural changes refer to the group of changes that result from local demographic, social, and economic changes, thus not directly caused by energy-related policies. In planning terms, this would correspond to the "business as usual" or "reference" scenario, i.e., the potential evolution of the system between two distinct moments in time, considering that no policies would be implemented in this timeframe. This comprises demographic, social, and economic changes, as well as climate variability.

The full list of factors of change associated with the local socio-economic and cultural changes that have been identified is presented in Table 5, along with the respective characterizing indicators.


**Table 5.** Factors of change associated with local socio-economic and cultural changes with corresponding indicators.

#### 3.2.1. Demographic Evolution

Demographic changes include the population growth (due to migration and natural growth) and also changes in the age pyramid and in the share of population living in urban and rural areas. All these variations imply changes in the local energy system—in the overall energy use as well as in the structure of this use. For example, energy use in urban areas differ from that in rural areas, due to the differences in heating options available (natural gas vs. traditional wood) and the different household conditions (housing type and size). A summarized description of each of the factors of change associated with demographic changes is presented below.

Population growth refers to the changes regarding the size of local population that occur within the studied timeframe. This is assessed through the number of local inhabitants and their evolution over time.

Migration to urban areas translates the level of urbanization of a municipality and its evolution over time. A change in the level of urbanization has an impact in the local energy system, in terms of preferred transportation modes, average distances traveled, and even in the type of housing. The changes in the level of urbanization can be assessed through the share of local population living in rural/urban areas.

Age pyramid refers to the distribution of the population according to their age. A change in this distribution implies changes in the local energy system, caused by the distinct energy service needs. The differences are more significant between active and inactive citizens. Thus, it is assumed that this factor of change can be assessed by the quantification of the share of active population.

#### 3.2.2. Social Evolution

Social changes refer to the lifestyle changes of local inhabitants that occur naturally, not promoted by energy-related policies. These may include e.g., the increase of appliances ownership per capita and changes in the number of persons per household, in the residential sector, and a change in the average distance between home and work or school, in the transportation sector. Each of the factors of change associated with social evolution are described below.

Household conditions refer to lifestyle changes associated with the residential building. This includes changes in the average number of people per household (usually also linked with the evolution in families' typology) and changes in the average size of the household. These changes can imply significant changes in the energy use for space heating and cooling, lighting, and even cooking. Moreover, changes in appliances ownership must also be considered.

Commuting needs reflect the distance between home and work/school for the local inhabitants. A variation in this distance implies also a change in the daily transportation needs, and often even a change in the transportation modes that are used for commuting. It is assumed that this factor can be assessed through the average distance traveled per person per day, for commuting purposes.

Transportation habits refer to changes in the preferences for daily traveling, for both commuting and leisure purposes. These characteristics may change over time for multiple reasons—economic, social or other. It is assumed that a large share of these changes can be estimated based on the changes in car ownership of local population.

#### 3.2.3. Economic Evolution

Economic evolution refers to the economic situation of local inhabitants and also to the evolution of local economic activities. The effect of these changes on local energy system may be affected, even if not in a very significant way, by local context. For instance, the increase in average income of local population leads to an increase in energy use; however, the growth rate may depend on the level of awareness of the local inhabitants towards climate change issues. The factors of change associated with the economic evolution at the local level are described below.

Overall growth of local economy refers to changes in the level of activity, i.e., the overall value of the local economy. Changes in the local gross domestic products (GDP) are a good indicator for the variation in the level of activity.

Structure of local economy refers to the weight that different economic sectors (or subsectors) represent to the overall local economy. Describing the structure of economy of a specific municipality would consist on identifying the sectors responsible for a larger share of gross value added or for employing more workers and which sectors are less relevant.

Economic situation of households translates the financial situation of individuals. This includes the employment situation, which has a significant impact in the local energy system given the difference in the daily habits between employed and unemployed citizens, and the individuals purchase power. On the latter, there are several studies that demonstrate the relation between private purchase power and energy use.

#### 3.2.4. Climate Variability

Finally, the effects of differences in climate conditions have also been included in the natural evolution of the local energy system. These vary from year to year, and temperature differences lead to changes in energy needs (mostly for heating and cooling), which must be considered. It is assumed that these changes can be assessed through the variations in annual heating and cooling degree days.

#### *3.3. Higher-Level Governance Framework*

Higher-level governance framework refers to the characteristics of international, national, and regional governance frameworks existing within the studied timeframe and focusing on energy and climate issues. Considering the variables identified in the review, it is assumed to be appropriate to separate technology evolution from specific energy and climate policies.

The full list of factors that characterize higher-level governance framework is presented in Table 6, along with the respective characterizing indicators.


**Table 6.** Factors of change related to higher-level governance framework with corresponding indicators.

#### 3.3.1. Technology Evolution

Technology evolution can be prompted by energy related policies (targeting energy efficiency and others) but it may also be market driven (following industry needs or an attempt to reduce costs). The evolution of the energy supply sector can be seen as a mix between natural technology evolution and an effect of energy and climate policies. The energy efficiency improvement on conversion plants is an example of the former, while the preference for renewable energy sources to produce electricity, due to fee-in tariffs in the technological maturation phase, is an example of the latter.

Technical evolution refers to the technology development that occurs naturally, i.e., excluding development that is promoted by energy efficiency policies and other related instruments. This evolution can be assessed by the variation in the energy intensity for each activity sector, at national or even international level.

Market evolution corresponds to the changes in the market that lead to changes in the technology choices by both individuals and companies. This can be translated into changes in technology costs, as well as the presence of market barriers that hamper the shift towards more efficient solutions. The variation in technology costs can be characterized by, e.g., the over cost of A++ appliances compared to a B equivalent. The presence of market barriers and respective evolution over time can be assessed by the quantification of the energy efficiency gap in different moments in time.

Energy supply sector represents the changes in the national and regional energy transformation and transport activities. It includes the evolution in electricity and heat production and transport, as well as the changes in the refinery of fossil fuels. Thus, it comprises the changes in the fuel-mix of electricity and heat generation, i.e., the increase/decrease in the share of renewable energy sources that are used in electricity production, as well as the changes in the efficiency of the energy supply sector (electricity and heat production, but also oil products transformation).

#### 3.3.2. National and International Energy-Related Policies

On what concerns the factors of change associated to energy and climate policies, two distinct groups of factors were found: the ones that characterize the existing energy and climate related policies, e.g., in terms of implementation and specific targets, and the ones that characterize the existing institutional framework regarding the promotion and support of local actions towards climate change mitigation. It is also important to mention the influence that the local context may have on the overall effect of these factors of change on the local energy system. The effects of national and international policies depend on the reaction of targeted actors to the policy instruments put in place. In turn, their reaction depends on their level of education and awareness as well as cultural habits and financial situation. All the identified factors associated with energy and climate policies are briefly described below.

Existing policies correspond to the policies that are implemented at regional, national, and international level that have a climate and/or energy purpose. These may impact positively or negatively local GHG emissions, depending on the coherence between local policies and the ones defined at higher levels of governance. Here, it is assumed that the public and private investments on environmental management and protection could be a good indicator of the scale of existing policies at governance levels other than local.

Energy costs refer to the changes in the average cost of the different energy products and the impact that these may have in the local energy system. Energy prices are usually established at the national and/or regional level, and thus are included within the changes associated with the higher-level governance framework. The average unitary price per energy product in the consumer sector, as well as the respective evolution over time, is considered a good indicator to characterize this factor of change.

Coherence with local policies refers to whether national and regional policy objectives are aligned with those that were defined at the local/municipal level. Here, the quantification of this alignment is not feasible, thus a yes or no question is suggested.

Support of local action refers to the existence of initiatives promoted at national and international level that promote actions taken at the local level. This may include technical support for the development of action plans, and energy dedicated actions, and economic incentives and support (as low-rate loans). A possible indicator to characterize the level of support could be the number of initiatives promoted at higher levels of governance that support local action.

#### *3.4. Local Climate Change Mitigation Actions*

Finally, local climate change mitigation actions include any policy actions that are promoted at the local level, by (or with the support of) local authorities, with the aim of reducing local GHG emissions. These include actions taken under a formal and wider plan and also isolated actions with a very specific target (e.g., the improvement of energy efficiency in public lighting). Here, the factors of change that have been identified are mostly related with the characterization of these actions, in terms of objectives and targeted actors as well as in what concerns implementation procedures. The factors that were identified as important for the characterization of local climate change mitigation actions are briefly described below and presented in Table 7.


**Table 7.** Factors of change associated with local climate change mitigation actions with corresponding indicators.

Actors' involvement refers to the group of actors that were/are involved in the local climate change mitigation actions. A distinction is made between the actors that are involved in the planning process (where actions are defined) and the ones involved in the implementation stage. The involvement of different stakeholders in the planning and implementation of climate and energy actions is considered to be beneficial for the respective acceptance and adherence by local citizens. The characterization of this involvement may be achieved by the identification of the groups of actors that were involved (e.g., local authority, local/regional energy agency, external consultant, local NGOs, citizens, etc.) and their level of involvement (co-creation, draft version public consultation, etc.).

Timeframe refers to the time span of the actions being planned/implemented. Throughout the review, the consideration of the time span of the ex-ante planning exercise (when existing) was also noted as relevant for the assessment of the outcomes of local actions. This is translated into the beginning and end years of the actual (or projected) action implementation.

Targets refer to the objectives that were defined at the local level. These include the level of GHG emissions reductions that is aimed, the wider vision concerning the local energy system, and whether there are other goals besides GHG emissions (such as creation of jobs or decrease in energy costs).

Current implementation status refers to the degree of implementation of the wider action plan (if existing) and of the different individual actions that are projected. It can be assessed as a percentage of the overall implementation target. For instance, considering an action that implied the installation of solar thermal systems for water heating in 100 residential buildings, the degree of implementation would correspond to the identification of the share of systems that were already installed.

Estimated impact corresponds to the projected outcome of the proposed actions, in terms of both energy use and GHG emissions. When a mid-term assessment is also performed, the quantification of the already achieved reductions could also be performed.

Budget refers to the level of investment needed to implement the proposed actions. Here, it is important to distinguish between the overall projected budget and the investments made thus far.

Financing sources correspond to the entities responsible for financing the actions' implementation. This may comprise different entities, both public and private. The characterization of the financing sources could be performed by the identification of the type of entities financing the local actions (local authority's own resources, national funds and programs, EU funds and programs, other), and the level of investment of each of the different entities.

Actions characterization corresponds to the characterization of the individual actions planned and/or implemented at the local level, in terms of: (1) area of intervention; (2) type of technical measure; and (3) policy instrument. Area of intervention is the sector and subsector of the local energy system which will be impacted by the action (residential, services, transports, etc.). Then, the technical measure refers to the type of change that is aimed at, as well as the degree of change that is intended. For instance, the shift of 10% of private cars passenger-kilometers to soft modes of transportation could be seen as a technical measure. Finally, the choice of policy instrument corresponds to the type of mechanism that is used to achieve the desired changes in the system. Using the same example, this could comprise the creation of dedicated bike lanes and/or the implementation of a congestion charge. When assessing the policy instruments, it is also useful to identify the level of financial support that is provided to the adopters of the measures (i.e., if there are economic incentives to change) and which actors are involved in their implementation.

Similar to what was said regarding the effects of higher-level policies, the influence of local context on the contribution of local actions cannot be disregarded. Local inhabitants' characteristics may influence the outcomes of policy actions in a significant manner. Moreover, local authorities' characteristics are an important feature on the definition and implementation of local actions. Thus, even if local actions are similar in their characteristics, their outcome can be very different if the context in which they are implemented is not the same.

#### **4. Discussion**

There are several studies that analyze the changes in local GHG emissions and potential explanatory factors. However, it was noted throughout the review that existing work is usually constrained to the consideration of a specific set of factors of change, disregarding others that may also have a significant impact on local emissions. For instance, it is common to narrow the focus to the actions that are being taken at the local level, disregarding the influence of higher-level governance framework in which they are being implemented and that may have confounding effects.

Such non-consideration of other potential factors of change hinders the accuracy of the results, as the effects calculated and associated with the considered factors may turn out to be either under- or overestimated. Indeed, the evaluation of the effects of local climate change mitigation actions will not lead to accurate results unless the remaining factors of change are taken into account. For example, in [5], the impact of local climate actions was inconclusive, partly due to the non-consideration of factors associated with the changes in local socio-economic and cultural conditions. Moreover, a better understanding of all potential factors of change of local GHG emissions can lead to better planning at local as well as regional and national levels. The interaction between the policies implemented at different governance levels could be improved by a comprehensive consideration of all factors of change and their interactions.

It is worth noting the overly simplified consideration of the local context in most of the existing work. Some of the reviewed documents refer to its importance, and how the local context specificities can be determinant for the success of local actions implementation, yet quantification methods tend to be overlooked. The systematic identification of the factors of change associated with local context may thus be a step forward for their consideration in future empirical ex-post evaluation studies and/or ex-ante planning exercises.

#### **5. Conclusions**

This paper is dedicated to a comprehensive and systematic identification and characterization of the factors that are relevant for the assessment of the observed changes in GHG emissions associated with local energy systems. This is an important contribution for the assessment of the actual effects of local climate change mitigation actions, by establishing the bases for developing analytical frameworks.

The identification was performed through a systematic literature review process. The complete list of variables encountered in the review was reduced using eligibility criteria, as well as three screening filters (explicit repetition, implicit repetition, and dependent variables). A list of 107 variables that were grouped into 40 factors of change was finally obtained.

The factors of change were grouped into four subsets, corresponding to three sources of change plus the local context. The sources of changes subsets are: local socio-economic and cultural changes; higher-level governance framework; and local climate change mitigation actions. The local context corresponds to a group of variables that, despite not leading directly to changes in local GHG emissions, have influence on how the remaining factors impact the local energy system—e.g., local inhabitants' level of awareness and commitment towards environmental problems or local authorities' legal scope of action.

Even if the identification of the major categories has been previously accomplished in academic studies, the major advancement of this work corresponds to their disaggregation into individual factors of change and respective characterization.

The factors identified correspond to all the factors that may potentially lead to changes in local GHG emissions. The relevance of the effects associated with each factor of change will depend on the local system in question, as well as on the timeframe addressed. Moreover, even acknowledging the difficulty of empirical validation (and quantification) of each factor individually, the obtained results comprise an advancement from what was possible thus far, providing a basis for more accurate models and analyses of the evolution of the energy systems at the local level.

These results can be used to build an analytical framework to characterize and model the evolution of a local energy system, assessing the individual effect of each of the different explanatory factors. This will allow for more accurate energy planning at the local level, taking into account the potential interactions with the external factors of change (i.e., the factors that are not controlled by local authorities). For instance, it may be used for ex-post evaluation of local climate change mitigation actions, allowing for the disaggregation of the observed changes in local GHG emissions into the causes and the isolation of the effects that were solely associated to local policies and actions.

The outcome of this work can be seen as a fundamental step towards a comprehensive assessment of the evolution of local energy systems. Such assessment is crucial for a better evaluation of the actual contribution of local climate-related policies and actions towards climate change mitigation.

**Author Contributions:** Conceptualization, I.A. and V.L.; Formal analysis, I.A. and V.L.; Funding acquisition, I.A.; Investigation, I.A.; Methodology, I.A.; Supervision, V.L.; Validation, I.A. and V.L.; Visualization, I.A.; Writing—original draft, I.A.; and Writing—review and editing, I.A. and V.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Fundação para a Ciência e Tecnologia (FCT—Portugal), grant number PD/BD/105862/2014, in the frame of the MIT Portugal Program. INEGI (LAETA—Grupo de Energia e Ambiente) provided financial support with the article processing charges.

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

#### **Appendix A**

This appendix presents in detail the process of identification and screening of the factors of change that are relevant for local GHG emissions. Table A1 presents the full list of variables identified in the literature. The four following tables are dedicated to the screening process: Table A2 presents the variables that were excluded for non-applicability; Tables A3 and A4 present the variables that were eliminated due to explicit and implicit repetition, respectively; and Table A5 presents the variables that were excluded to avoid dependent variables.

**Table A1.** Full list of variables identified in the systematic literature review process and respective sources.



**Table A1.** *Cont.*



**Document Code Variable** IAEA (2006) 12.040 Electricity penetration for appliances in urban households 12.041 Specific fossil fuel consumption per urban dwelling for lighting and non-electric appliances 12.042 Penetration of various energy forms into space heating in urban households 12.043 <sup>E</sup>fficiency of various fuels use, relative to that of electricity use, for space heating in urban households 12.044 Coefficient of performance of heat pumps for space heating in urban households 12.045 Penetration of various energy forms into water heating in urban households 12.046 <sup>E</sup>fficiency of various fuels use, relative to that of electricity use, for water heating in urban households 12.047 Coefficient of performance of heat pumps for water heating in urban households 12.048 Approximate share of water heating demand in urban households that can be met with solar installations 12.049 Penetration of various energy forms into cooking in urban households 12.050 <sup>E</sup>fficiency of various fuels use, relative to that of electricity, for cooking in urban households 12.051 Approximate share of cooking demand in urban households that can be met with solar installations 12.052 Share of air conditioning demand of urban households that can be met with electricity 12.053 Coefficient of performance of electric air conditioning in urban households 12.054 Coefficient of performance of non-electric air conditioning in urban households 12.055 Fraction of rural dwellings in areas where space heating is required 12.056 Degree days for rural dwellings 12.057 Fraction of rural dwellings per type 12.058 Average size of rural dwellings by type 12.059 Fraction of floor area that is actually heated in rural areas per dwelling type 12.060 Specific heat loss rate by rural dwelling type 12.061 Share of rural dwellings with air conditioning, by dwelling type 12.062 Specific cooling requirements by rural dwelling type 12.063 Specific energy consumption for cooking in rural dwellings 12.064 Share of rural dwellings with hot water facilities 12.065 Specific energy consumption per rural dwelling for electric appliances 12.066 Electricity penetration for appliances in rural households 12.067 Specific fossil fuel consumption per rural dwelling for lighting and non-electric appliances 12.068 Penetration of various energy forms into space heating in rural households 12.069 <sup>E</sup>fficiency of various fuels use, relative to that of electricity use, for space heating in rural households 12.070 Coefficient of performance of heat pumps for space heating in rural households 12.071 Penetration of various energy forms into water heating in rural households 12.072 <sup>E</sup>fficiency of various fuels use, relative to that of electricity use, for water heating in rural households 12.073 Coefficient of performance of heat pumps for water heating in rural households 12.074 Approximate share of water heating demand in rural households that can be met with solar installations IAEA (2006) 12.075 Penetration of various energy forms into cooking in rural households 12.076 <sup>E</sup>fficiency of various fuels use, relative to that of electricity, for cooking in rural households 12.077 Approximate share of cooking demand in rural households that can be met with solar installations

**Document Code Variable** IAEA (2006) 12.078 Share of air conditioning demand of rural households that can be met with electricity 12.079 Coefficient of performance of electric air conditioning in rural households 12.080 Coefficient of performance of non-electric air conditioning in rural households 12.081 Share of service sector in the total active labor force 12.082 Average floor area per employee in the Service sector 12.083 Share of Service sector floor area requiring space heating 12.084 Share of Service sector floor area requiring space heating that is actually heated 12.085 Specific space heat requirements of Service sector floor area 12.086 Share of air-conditioned Service sector floor area 12.087 Specific cooling requirements in the Service sector 12.088 Energy intensity of motor fuel use per Services' subsector 12.089 Energy intensity of electricity specific uses per Services' subsector 12.090 Energy intensity of thermal uses (except space heating) per Services' subsector 12.091 Penetration of various energy forms into space heating in Service sector 12.092 Penetration of various energy forms into other thermal uses in the Service sector 12.093 <sup>E</sup>fficiency of various fuels use, relative to that of electricity use, for thermal uses in Service sector 12.094 Coefficient of performance of heat pumps in space heating in Service sector 12.095 Share of low-rise buildings in the total Service sector floor area 12.096 Approximate share of thermal uses in the Service sector that can be met by solar installations 12.097 Share of air-conditioning that can be met with electricity 12.098 Coefficient of performance of electric air conditioning in Service sector 12.099 Coefficient of performance of non-electric air conditioning in Service sector IEA (2014) 13.001 Population 13.002 Floor area per capita 13.003 Persons per household 13.004 Appliances ownership per capita 13.005 Heat per floor area 13.006 Domestic how water energy per capita 13.007 Cooking energy per capita 13.008 Electricity for lighting per floor area 13.009 Energy per appliance 13.010 Passenger-kilometers 13.011 Share of total passenger-kilometers by mode 13.012 Energy per passenger-kilometer by mode 13.013 Tons-kilometers 13.014 Share of total ton-kilometers by mode 13.015 Energy per ton-kilometer by mode 13.016 Services value-added 13.017 Energy per value-added (services) 13.018 Industry subsectors value-added 13.019 Share of total value-added by subsector 13.020 Energy per value-added by subsector 13.021 Agriculture value-Added 13.022 Share of total value-added by subsector 13.023 Energy per value-added by subsector 13.024 Changes in CO2 intensity 13.025 Changes in input coefficients 13.026 Changes in the composition of final demand 13.027 Changes in the level of final demand (economic growth)

**Table A1.** *Cont.*


**Table A1.** *Cont.*



**Table A1.** *Cont.*


**Table A1.** *Cont.*


**Table A2.** Variables excluded for non-applicability in the systematic literature review process.



**Table A2.** *Cont.*


**Table A3.** Variables excluded for explicit repetition in the systematic literature review process.



**Table A3.** *Cont.*

**Table A4.** Variables excluded for implicit repetition in the systematic literature review process.





**Table A5.** Variables excluded to avoid the presence of dependent variables in the systematic literature review process.




#### **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* **Wind Power Cogeneration to Reduce Peak Electricity Demand in Mexican States Along the Gulf of Mexico**

**Quetzalcoatl Hernandez-Escobedo 1, Javier Garrido 1, Fernando Rueda-Martinez 1, Gerardo Alcalá <sup>2</sup> and Alberto-Jesus Perea-Moreno 3,\***


Received: 14 May 2019; Accepted: 10 June 2019; Published: 18 June 2019

**Abstract:** The Energetic Transition Law in Mexico has established that in the next years, the country has to produce at least 35% of its energy from clean sources in 2024. Based on this, a proposal in this study is the cogeneration between the principal thermal power plants along the Mexican states of the Gulf of Mexico with modeled wind farms near to these thermal plants with the objective to reduce peak electricity demand. These microscale models were done with hourly MERRA-2 data that included wind speed, wind direction, temperature, and atmospheric pressure with records from 1980–2018 and taking into account roughness, orography, and climatology of the site. Wind speed daily profile for each model was compared to electricity demand trajectory, and it was seen that wind speed has a peak at the same time. The amount of power delivered to the electric grid with this cogeneration in Rio Bravo and Altamira (Northeast region) is 2657.02 MW and for Tuxpan and Dos Bocas from the Eastern region is 3196.18 MW. This implies a reduction at the peak demand. In the Northeast region, the power demand at the peak is 8000 MW, and for Eastern region 7200 MW. If wind farms and thermal power plants work at the same time in Northeast and Eastern regions, the amount of power delivered by other sources of energy at this moment will be 5342.98 MW and 4003.82 MW, respectively.

**Keywords:** wind farm; thermal power plants; peak electricity demand; Gulf of Mexico

#### **1. Introduction**

One of the world's biggest problems is that population always is increasing and some countries have over dependence on energy generation. In the world, primary energy consumption growth averaged 2.2% in 2017, up from 1.2% last year and the fastest since 2013. This compares with the 10-year average of 1.7% per year [1]. In 2016, generation from combustible fuels accounted for 67.3% of total world gross electricity production (of which: 65.1% from fossil fuels; 2.3% from biofuels and waste), hydroelectric plants (including pumped storage) provided 16.6%; nuclear plants 10.4%; geothermal, solar, wind, tide, and other sources 5.6%; and biofuels and waste made up the remaining 2.3% [2]. It is known that energy supply side and demand side is essential to the correct usage of energy [3], and this is important to nations because it depends on the effectiveness of its power systems and the availability of power supply at the load centers [4]. It is important to define the concepts of load and electricity demand, on the one hand to understanding load of an electric power distribution as the final stage of the system that converts electric energy in another form of energy; respecting electricity demand is known as the electric power related on a period of time, that uses load to work.

In Mexico, electricity demand is considered decisive in social and economic development and consequently in the improvement of the economic conditions. The maximum demand for electrical energy occurs mainly during the summer season, between May and September, when the heat is more severe in the north and southeast of the country. As a result, the use of ventilation equipment is increased to maintain the temperature at more comfortable levels and the values necessary for the proper functioning of certain appliances are regulated. In 2017, the energy independence index, which shows the relationship between production and national energy consumption, was equivalent to 0.76. This result implies that the amount of energy produced in the country was 24.0% less than that which was made available in the various consumer activities in the national territory. In the course of the last ten years, this indicator decreased in an average annual rate of 5.0% [5]. An alternative will be the use of renewable energy, the Transition Energetic Law establish a renewable energy minimal participation in electricity generation of 25%, 30% and 35% in 2018, 2021, and 2024, respectively [6]. In the energetic mix, renewable energies have a participation of 14%, the technologies most used are hydraulic, and wind energy contributed with 3% or 38.23 PJ in 2017 [7].

Peak electricity demand is a global policy concern which creates transmission constraints and congestion, and raises the cost of electricity for all end-users [8]. Also, a considerable investment is required to upgrade electricity distribution and transmission infrastructure, and build generation plants to provide power during peak demand periods [9]. This is important because usually service suppliers charge a higher price for services at peak-time than for off-peak time to compensate for the costly electricity generation at peak hours [10]. Thus, if peak demand is reduced, electric system will be benefited. In addition to these benefits a study done by the authors of Reference [11] determine that it can eliminate the need to install expensive extra generation capacity such as combustion turbines for peak hours which are less than a hundred hours a year. With this information regions will be considered where demand complies with the conditions to have the necessary hours. Another benefit is for consumers, because if they know the period of time of peak demand and decrease their consumption, they can avoid higher electricity prices. An option to reduce peak demand is installing peak plants but the economic benefit does not justify the investment. Usually these peak plants use fossil fuel and the result is more generation of greenhouse gases [12].

The use of wind energy to support the electric demand mainly during its peak period is an alternative to reduce fossil fuel combustion and greenhouse gas emissions. Wind is intermittent and uncertain, which represents a problem in power systems [13]. A study done by the authors of Reference [14] established that due to the unpredictable nature of wind energy and non-coincidence between wind units output power and demand peak load, wind units are deemed as an unreliable source of energy. In this study a stochastic mathematical model was developed for the optimal allocation of energy storage units in active distribution networks in order to reduce wind power spillage and load curtailment while managing congestion and voltage deviation.

In New York State the intraregional effects were illustrated by quantifying the net load, net load ramping, operating reserve, and regulation requirements, the study found out that only at wind capacities exceeding 100% of the average statewide load does the wind-generated electricity meet significant portions of the distant demands [15].

An improved energy hub system combined with cooling, heating, and power integrated with photovoltaic and wind turbine, demonstrated that operation cost and carbon emission were reduced by 3.9% and 2.26%, respectively [16]. In Canada a battery storage system to minimize the energy drawn from the public grid was proposed This storage system was charged from wind turbines and the energy stores were discharged to a park when the wind park power dropped below 0 kW and the storage system was able to offset 17.2 MWh but the financial gain was insufficient to offset the net energy losses in the storage system [17].

About cogeneration, in China Reference [18] established the relationship between a wind farm and a thermal power plant, in this study the difference between wind power profits and thermal power is only analyzed in terms of cost differences, and no comparison of other data is involved. According to this, the specific example of Gansu Province not only shows the difference in earnings before and after the peak-to-peak adjustment of wind power, but also indicates that the large-scale auxiliary adjustment of wind power can ensure that thermal power generators can participate in deep-peaking with reasonable returns; Reference [19] analyses and compares several existing kinds of distributed energy storage for improving wind power integration. Then considering a mass of cogeneration units in north China, control strategy for wind power integration is proposed to reduce the peak regulation capacity for wind power integration. As a result, the peak–valley load difference can be equivalently reduced to relieve peaking pressure for wind power integration.

In Northern China a proposal were made suggesting that if the energy carrier for part of the end users space heating is switched from heating water to electricity (e.g., electric heat pumps can provide space heating in the domestic sector), the ratio of electricity to heating water load should be adjusted to optimize the power dispatch between cogeneration units and wind turbines, resulting in fuel conservation, they found out that with the existing infrastructures are made full use of, and no additional ones are required [20]. Bexten et al. [21] investigate a system configuration, which incorporates a heat-driven industrial gas turbine interacting with a wind farm providing volatile renewable power generation. The study quantifies the impact of selected system design parameters on the quality of local wind power system integration that can be achieved with a specific set of parameters. Results show that the investigated system configuration has the ability to significantly increase the level of local wind power integration.

In this study a methodology is proposed, which consists in locating thermal power plants along the five states of the Gulf of Mexico (Tamaulipas, Veracruz, Tabasco, Campeche, and Yucatan) and evaluating wind resources next to these plants with the objective to determine the power output generated and their contribution to reduce the peak electricity demand. For the geographic position of Mexican states along the Gulf of Mexico, see Figure 1.

**Figure 1.** Mexican states along the Gulf of Mexico.

#### **2. Material and Methods**

#### *2.1. Data*

Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) provides re-analysis data such as: wind speed and wind direction at 50 m height, temperature at 2 m and 10 m, and atmospheric pressure with records from 1980–2018. The resolution of every point downloaded is 50 km covering all the Mexican states along the Gulf of Mexico including offshore points, which can all be observed in Figure 2. Several studies have used MERRA-2 to determine wind characteristics, as did Reference [22] in the Central Californian Coast, assessing offshore conditions. A study about wind and rainfall areas of tropical cyclones making landfall over South Korea was examined for the period 1998–2013 using MERRA-2 data. It was determined that composite analyses of the cases of strong and weak vertical wind shear confirm that the increase of rainfall area is related to the asymmetric convection (rising/sinking motion in the downshear-left/upshear-right side) induced by the vertical wind shear [23]. An investigation on prediction of wind power because of climate change was made using MERRA-2 among other re-analysis data [24]. As can be seen, MERRA-2 data have been used with excellent results on wind assessment.

**Figure 2.** Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) data position.

In total, 142 MERRA-2 points, 61 in Tamaulipas, 31 in Veracruz, 9 in Tabasco, 24 in Campeche and 17 in Yucatan.

#### *2.2. Thermal Power Plants*

There are 22 thermal power plants in Mexico that generate electricity. These plants work on a combined cycle with gas turbines. This generation contributes 55.6% of all electricity in the country [25]. Along the Gulf of Mexico, the most important plants are located as follows: two in the state of Veracruz, one called Adolfo Lopez Mateos with capacity of 2263 MW and Dos Bocas with capacity 350 MW. In the state of Tamaulipas there are two, Rio Bravo and Altamira, with capacity of 800 MW and 1039 MW respectively [26]. Figure 3 shows the geographic position of these thermal power plants. Unfortunately, there are no thermal plants located in the states of Tabasco, Campeche, and Yucatan.

**Figure 3.** Thermal power plants.

#### *2.3. Wind Assessment*

Wind characterization is essential to determine its power at a given site. Its assessment depends mainly on the orography, roughness, and generalized climatology.

Statistically it can be characterized as a probability density function (PDF) that can be interpreted as the probability that the random variable *x* lies in a differential range, *dx*, about a value *x* is *f(x)dx*. More specific statements about the probability that the unsteady variable *x*, lies in a particular range *a* ≤ *x* ≥ *b*, this expression is given by Equation (1).

$$P = \left(a \le x \le b\right) = \int\_{a}^{b} f(x)dx,\tag{1}$$

Another PDF used in wind assessment is the Weibull distribution for the wind speed (*u*) that has two parameters: shape parameter (*k*) and scale parameter (*c*). This is expressed in Equation (2):

$$f(u) = \frac{k}{c} \left(\frac{u}{c}\right)^{k-1} e^{-\left(\frac{u}{c}\right)^k},\tag{2}$$

The power in the wind is the product of mass flow rate through turbine blades ρ*UA*, and the kinetic energy per unit mass in the wind, *U2*/2. The average power in the incoming wind is presented in Equation (3) [27].

$$P = \frac{1}{2} \rho A \mathcal{U}^3 \mathcal{C}\_{P'} \tag{3}$$

where ρ is the density of air, *CP* is the power coefficient, *A* is the rotor swept area, *U* is the wind speed and *P* is the wind turbine power output.

#### *2.4. Peak Electrical Demand*

In Mexico the National Energy Control Center (NECC) divided the country in seven regions. This division is called National Interconnected System [28], and is based in two electric regions: Northeast and Eastern [29]. These regions can be seen in Figure 4.

**Figure 4.** Regions analyzed.

NECC shows information about demand forecasting, net demand, and total demand. It uses three methodologies to calculate demand forecasting: moving average that is optimal for random or leveled demand patterns where the impact of historical irregular elements is to be eliminated through an approach of periods of up to 7 days previously. This methodology is the most used by NECC and is given by Equation (4)

$$\hat{X}\_t = \frac{\sum\_{i=1}^n x\_{t-1}}{n},\tag{4}$$

where *<sup>X</sup>*ˆ*<sup>t</sup>* is the energy demand average in period *<sup>t</sup>*; *Xt*−<sup>1</sup> is the real demand for the periods prior to *<sup>t</sup>* and the number of observations *n*.

The second method is weighted moving average, a variation of moving average. While in moving average all data has the same weight, in weighted moving average each data is assigned with a specific weight, always considering that total sum will be 100% and is expressed as Equation (5):

$$\hat{X}\_t = \sum\_{i=1}^n c\_i x\_{t-1} \,\prime \tag{5}$$

where *<sup>X</sup>*ˆ*<sup>t</sup>* is the energy demand average in period *<sup>t</sup>*; *ci* weighted factor; *xt*−<sup>1</sup> is the real demand for the periods prior to *t* and the number of observations *n*.

The third method is called multiple linear regression, where both independent and dependent variables working with the independent variable variation to forecasting the dependent variable. Multiple linear regression is given by Equation (6):

$$D\_{\mathbf{F}} = \beta\_0 + \beta\_1 \mathbf{x}\_1 + \dots + \beta\_n \mathbf{x}\_n + \varepsilon,\tag{6}$$

where *DF* is Demand Forecasted; *x*1, ... *x*<sup>n</sup> independent variables; β*<sup>0</sup>* ... β*<sup>n</sup>* coefficients which are calculated by least squares; *n* number of data and ε is a random error term.

Several studies about peak demand where demand consumers participate have established two typical types of programs. The first one is an incentive-based program the end users are encouraged by the incentive payments to reduce their load demands, and thus to help improve the feasibility and stability of the power system [30]. This program usually includes direct load control and demand bidding [31]. The second one is price-based, and in this case final consumers are encouraged to configure their energy usage patterns according to electricity demand curve. These price-based programs consist in observing critical peaking pricing, time of use pricing, and time pricing [32,33].

Figure 5 presents two demand curves, each one by region. These demand curves have been developed with the aim to estimate the current hourly electricity demand profiles with a temporal resolution of 1 h based in 2018.

(**b**)

**Figure 5.** (**a**) Northeast electricity demand and (**b**) Eastern electricity demand.

As seen in Figure 5, the behavior of electric demand varies among regions, Figure 5a at Northeast the peak appears in range of 16:00–18:00 h and has another peak during 22:00–24:00 h. This is because the consumption increases in the first peak due to domestic use and the second electricity peak is due to industrial sector and in Figure 5b the Eastern region presents a peak demand during 21:00–23:00 h, mainly for industrial activities.

These demand curves will be used and compared to wind daily profile to consider it as wind power.

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

Along the Gulf of Mexico 142 MERRA-2 points have been assessed with wind speed information. 4 sites around thermal power plants have been chosen and the amount of power output generated has been calculated to reduce peak electricity demand.

#### *3.1. Wind Resource Assessment*

Wind has been converted as resource along the Gulf of Mexico; this means that mainly both roughness and orography have been taken into account.

#### 3.1.1. Rio Bravo, Altamira, Tuxpan, and Dos Bocas

Four MERRA-2 points were located near thermal power plants, their geographic information is presented in Table 1.

**Table 1.** Geographic information for Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) and thermal power plants.


In Appendix A Tables A1–A4 hourly mean wind speed are presented by month for MERRA-2 points. In each one of them wind speed behavior can be appreciated.

A typical wind year can be compared to electricity demand to identify its trajectory, so that it could be seen if in the peak demand there existed a peak wind. To probe it, Figures 6 and 7 represent both electricity demand and wind speed of Rio Bravo, Altamira, Tuxpan, and Dos Bocas respectively.

As can be observed in Figure 6, wind speed trajectory has an increase in the last hours as the electricity demand does as follows: Figure 6a 6.8 m/s to 7.6 m/s and Figure 6b5m/s to 6.5 m/s both during 21:00 h to 24:00 h. Although Figure 6b does not follow the entire trajectory, at the end of the day the last 4 h (21:00–24:00) when the demand begins to grow, wind also does. In Figure 7a 4.5 m/s to 6.5 m/s and Figure 7b5m/s to 6 m/s both between 20:00 h to 24:00 h. This could be considered to model a potential source of power.

Determining wind speed features correctly means that the model that will be designed will have more accuracy. In Table 2, statistical information at the points studied is presented.

**Figure 6.** (**a**) Northeast demand compared to wind speed at Rio Bravo and (**b**) Northeast demand compared to wind speed at Altamira.

**Figure 7.** (**a**) Eastern electricity demand compared to wind speed at Tuxpan and (**b**) Eastern electricity demand compared to wind speed at Dos Bocas.


**Table 2.** Wind resource statistics.

#### 3.1.2. Microscale Model

One of the critical criteria to place wind turbines is to assess wind resource and determine if the site have features as transmission lines, substations, and roads, among others.

A microscale model has been done on WAsP software including roughness and orography features. In Figures 8 and 9 the four sites are shown microscale modeled. Each one has an array of seven wind turbines placed to avoid wake effect. Wind turbines generators (WTG) are selected according to wind features at the site.

Modeled areas show different types of wind resource, and the resolution is 10 km. The areas of high wind energy potential are better identified. An advantage of microscale modelling is to extrapolate onsite measurements to the prospective turbine sites to modify the background atmospheric conditions and account for topographic induced speedups, turbulence, and wind shear [34]. The *x* and *y* axes are represented in UTM scale, in Figure 8a the zone presents a wind power density between 261 W/m2 to 286 W/m2, WTG type GEWind 77 were placed to avoid wake effects; Figure 8b presents a zone with 170 W/m2, level curves show hills in this area and were harnessed to place WTG type BONUS 600 kW along it. In Figure 9a, this zone is in front of the sea, an area with 142 W/m2, wind turbines VESTAS V82-1650 kW were used; Figure 9b. In Dos Bocas there is a wind power density of 144 W/m2. Seven WTC type BONUS 600 kW were placed at the top of the hills.

(**a**)

**Figure 8.** (**a**) Rio Bravo and (**b**) Altamira.

(**a**)

**Figure 9.** *Cont*.

(**b**)

**Figure 9.** (**a**) Tuxpan and (**b**) Dos Bocas.

3.1.3. Power Output Generation and Electricity Peak Demand

In the northeast region, peak demand occurs twice per day, the first period at 17:00 h and the second one at 24:00 h. These peaks represent a higher price of electricity by consumers. Wind speed at Rio Bravo behaves similarly to electricity demand curve, wind daily profile has two peaks, and these at the same time. Meanwhile, wind speed profile in Altamira only has one peak at the end of the day, even so power can be extracted from the wind. The power contribution by two wind farms placed in Rio Bravo and Altamira in the state of Tamaulipas, respectively, is shown in Table 3.


**Table 3.** Wind farms power generation along the state of Tamaulipas.

In the state of Veracruz located at the eastern electricity region, wind has been analyzed near two thermal power plants. In this region, electricity demand has a peak that occurs at 22:00 h in both Tuxpan and Dos Bocas, and wind speed trajectory shows that it can contribute delivering electricity to the grid at the same time that peak demand occurs. Table 4 presents the power output generated of two wind farms at the state of Veracruz.

**Table 4.** Wind farms power generation along the state of Veracruz.


If the period of time and the wind speed are taken into account when peak demand is increasing, it will be able to calculate how much energy could be extracted in each point, see Table 5.


**Table 5.** Total power generation at peak demand.

As shown in Table 5, electric power generated by each wind farm could contribute to reduce peak demand. The electric power generated by region indicate an amount of 2657.02 MW for Northeast and 3196.18 MW for Eastern. This is the power delivered to the electric grid.

Wind resource assessment is fundamental to determine the site where wind can be used to generated electric power. Another variable that can be considered is the place, which is important taking into account electricity required, as mentioned in Reference [9]. In this work, wind farms were modeled near to thermal power plants and used their electric infrastructure. It was found the hours at day were both electricity demand and wind speed have their peaks, whit this information could benefit end-users as studied [12] that establish if consumers know the period of time were peak demand occurs could avoid higher electricity prices. This work presents a methodology to assess wind speed and calculate the power delivered to the grid, but is important to mention that if only wind speed is assessed and other variables such as roughness, orography, and climatology are not considered, wind cannot be considered as an unreliable source, as did Reference [15]. To avoid this problem, a proposal of cogeneration with thermal power plants has been made, in comparison with Reference [18] where a cost analysis was its main objective, we have analyzed wind data trajectory to reduce peak demand.

The proposal of connection to the power grid is based on the thermal power plant of Dos Bocas. The same connection is considered for the other thermal plants. In this case, the nominal power, that is, 350 MW, and the wind farm nominal power of 10.5 MW, is equivalent to 3% of the total substation capacity which indicates that the connection is possible. For Rio Bravo 800 MW, the connection will represent 1.3%, Altamira 1039 MW will be 1%, and Tuxpan 2263 MW. This connection represents 0.4% of the substation capacity.

Regarding the impact of new wind farms on power grids, Reference [35] established that power losses in distribution systems vary with numerous factors depending on the system configuration, such as level of losses through transmission and distribution lines, transformers, capacitors, and insulator, among others. In the case of this study, all the power losses were calculated in the original design. Appendix B shows a single line diagram where the transmission lines and the node modified with the connected wind farm can be seen.

#### **4. Conclusions**

The microscale model increased the accuracy of wind assessment to place wind farms. This is because when running the microscale model, several variables are taken into account such as roughness, orography, and climatology. In this study, four sites were assessed to model wind farms using MERRA-2 data, wind speed, wind direction, temperature, and atmospheric pressure, with records from 1980–2018.

The electricity that wind farms generate is delivered over transmission and distribution power lines. High-voltage transmission lines carry electricity over long distances to where consumers need it. More distance means more losses. The amount of power produced by each wind farm could contribute to the power produced by a thermal plant. In total for Northeast region they can deliver to the grid 2657.02 MW at peak demand. If there is 8000 MW demanded at this moment, the reduction will be at 5342.98 MW. At the Eastern region, the thermal power plant called Tuxpan is the biggest one in Mexico, so the amount of power in both wind farms and thermal plants is higher than Northeast. In this case, the total production is 3196.18 MW, and peak demand power delivered is 7200 MW if wind farms and thermal power plants work at the same time. The amount of power at this moment will be 4003.82 MW. There could be several techniques to contribute to decreased peak demand, energy storage technologies, and faster response times, but the most important is to nurture future generations.

The impact of the thermal power plants, once they have the contribution of wind energy, can be considered as follows: the efficiency for Rio Bravo and Altamira could be 31% and for Tuxpan and Dos Bocas 47%; the capacity factor, 36% and 56% respectively, however, the volatility and availability of wind in the entire year must be considered.

Mexico's main objective is to reach 35% of its electricity generated by clean sources in 2024. Employing this type of generation is important to contemplate both the trajectory of electricity demand and wind speed and to avoid a duck-curve that represents several renewable sources being contributed during the day but not at night, as solar or when wind does not has enough potential, e.g., Northeast curve has two peaks during the day, the first one peak is at 17:00 h and the second one is at 24:00 h, at this hour solar energy cannot be consider and if at this time does not flow the wind, the peak demand will increase and the prices will also do it.

**Author Contributions:** Q.H.-E. and A.-J.P.-M. have contributed to the theoretical approaches, simulation, experimental tests, and preparation of the article, and have contributed equally to this work. J.G., F.R.-M. and G.A. All authors have read and approved the final manuscript.

**Funding:** This research was funded by Ministry of Education Mexico, under the project Support to Research Groups, UV-CA-466.

**Acknowledgments:** Authors acknowledge to research group "Wind and solar energy" from Universidad del Istmo.

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

#### **Appendix A**


**Table A1.** Hourly mean wind speed in Rio Bravo Mexico.


**Table A2.** Hourly mean wind speed in Altamira Mexico.

**Table A3.** Hourly mean wind speed in Tuxpan Mexico.



**Table A4.** Hourly mean wind speed in Dos Bocas Mexico.

#### **Appendix B**

**Figure A1.** Single line diagram.

#### **References**


© 2019 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* **Seasonal Wind Energy Characterization in the Gulf of Mexico**

#### **Alberto-Jesus Perea-Moreno 1, Gerardo Alcalá <sup>2</sup> and Quetzalcoatl Hernandez-Escobedo 3,\***


Received: 21 November 2019; Accepted: 19 December 2019; Published: 23 December 2019

**Abstract:** In line with Mexico's interest in determining its wind resources, in this paper, 141 locations along the states of the Gulf of Mexico have been analyzed by calculating the main wind characteristics, such as the Weibull shape (*c*) and scale (*k*) parameters, and wind power density (WPD), by using re-analysis MERRA-2 (Modern-Era Retrospective Analysis for Research and Applications version 2) data with hourly records from 1980–2017 at a 50-m height. The analysis has been carried out using the R free software, whose its principal function is for statistical computing and graphics, to characterize the wind speed and determine its annual and seasonal (spring, summer, autumn, and winter) behavior for each state. As a result, the analysis determined two different wind seasons along the Gulf of Mexico;, it was found that in the states of Tamaulipas, Veracruz, and Tabasco wind season took place during autumn, winter, and spring, while for the states of Campeche and Yucatan, the only two states that shared its coast with the Caribbean Sea and the Gulf of Mexico, the wind season occurred only in winter and spring. In addition, it was found that by considering a seasonal analysis, more accurate information on wind characteristics could be generated; thus, by applying the Weibull distribution function, optimal zones for determining wind as a resource of energy can be established. Furthermore, a *k*-means algorithm was applied to the wind data, obtaining three clusters that can be seen by month; these results and using the Weibull parameter *c* allow for selecting the optimum wind turbine based on its power coefficient or efficiency.

**Keywords:** Weibull function; scale parameter; shape parameter; wind power density; seasons; Gulf of Mexico

#### **1. Introduction**

Global electricity demand grew 4% in 2018, almost twice than that for 2010, where renewables and nuclear power met most of the growth in demand [1]. According to the International Renewable Energy Association (IRENA) [2], 171 GW were added in 2018 worldwide after a strong growth in the last decade in renewable energy capacity. The total use of all renewables increased by 7.9%, where wind and solar energy contributed 84% of this total, and it is expected that in 2023, this increase for wind and solar energy will be 12.4% and 24%, respectively, in 2030, where the solar photovoltaic and wind power will be key energy sources since they are the energies with the highest growth, with the latter being one of the most profitable sources of energy in the world [3]. One of the most studied atmospheric parameters for decades is the direction of the wind [4]; nowadays this parameter is essential for the installation of a wind farm, where it is important that the wind turbines are not reducing their energy capacity due to poor design [5].

The role that Mexico will play on the renewable energy scenario is mainly outlined by its Energy Transition Law, which was recently reformed in 2015, establishing a minimum share of clean energy in the generation of electric power of 25% for 2018, 30% for 2021, and 35% for 2024 [6]. In 2017, there was a total electricity consumption of 260,051.895 GWh, out of which, renewables contributed around 77,907.2 GWh (30% of electric generation), where wind energy contributed 10,378 GWh (4% of electric generation), generated by 46 wind farms around Mexico [7].

Properly assessing wind resources for electricity generation implies knowledge of its behavior, which involves considering several variables that range from climate to the corrected wind turbine power curve (WTPC) selection [8]. Several studies in Mexico have located zones with suitable wind power, such as in the Baja California Peninsula, where it was found that the wind power density was above 400 W/m2 [9]. Another study investigated the social impact caused by the expansion of large-scale wind energy projects on the Isthmus of Tehuantepec [10]. Furthermore, in the state of San Luis Potosi, a study of the installation of a hybrid PV-wind power generation system for social interest houses was validated [11]. Also, a statistical methodology based on support vector regression for wind-speed forecasting located at La Ventosa, Oaxaca, Mexico, showed that forecasts made with their own method are more accurate for medium (5–23 h ahead) short-term wind speed forecasting (WSF) and wind-power forecasting (WPF) than those made with persistent and autoregressive models [12] and a wind power map of all of Mexico [13].

According to Mazzeo et al. [14], numerous studies have been done to classify patterns of weather variables, where some are related to cluster analysis (CA) [15,16] and other to the similarity of time-series wind vectors [17]. CA has been successfully applied to the regionalization of wind in complex terrain, such as those of the north of the Iberian Peninsula [18]. Furthermore, in the long-term, it will be useful to detect wind patterns through the analysis of their time series with fast Fourier transform (FFT) [19] or wavelet methods [20]. The use of the latter technique has proved useful for missing wind data arrangement for a long time series of wind data [21].

Therefore, the crucial step in wind assessment is to determine the wind resources, which depends on accurate wind-speed modelling [22]. One of these models is the probability density function (PDF) because it provides important wind speed distribution parameters and allows one to determine the Weibull parameters [23]. It is known that if wind speed follows the Weibull distribution with scale parameter (*c*), which has the same units as wind speed and the dimensionless shape parameter (*k*), the load and power density also follow the same distribution with shape parameters *k*/2 and *k*/3, respectively, and scale parameters *c*/2 and *c*/3, respectively [24]. To calculate the Weibull parameters, some processes have been developed, e.g., Saleh et al. [24] reviewed six kinds of numerical methods commonly used for estimating the Weibull parameters: the moment, empirical, graphical, maximum likelihood, modified maximum likelihood, and energy pattern factor methods. From the review, they found that if the wind speed distribution matched well with the Weibull function, the six methods were applicable; but if not, the maximum likelihood method performed best, followed by the modified maximum likelihood and moment methods. Furthermore, Aukitino et al. [25] assessed the wind speed and found out that moment method was the best fitting, Akda ˘g and Dinler [26] proposed a new method based on wind power density and mean wind speed called the power density (PD) method, Baseer et al. [27] estimated the Weibull parameters using a least-squares regression method (LSRM) and the Wind Atlas Analysis and Application Program (WAsP) algorithm, and Ozay and Celiktas [28] did an statistical analysis in Turkey using WAsP.

Wind power density (WPD) is one of the most important factors to consider when wind is assessed to generate power and depends on an accurate determination of the Weibull parameters. The literature shows the relationship between the WPD and the Weibull function. Katinas et al. [29] performed a study in Lithuania where the WPD was calculated after determining the Weibull parameters. In Spain, an evaluation of the WPD was done using the moment method for calculating the Weibull parameters [30]. Faghani et al. [31] extrapolated data at high altitudes, calculated the

Weibull parameters, and found out that variation of power density with time was significant; therefore, they divided the year in two periods, period I (spring and summer) and period II (autumn and winter).

The ideas of Faghani et al. [31] are considered for this study and are extended for seasonal analysis (spring, summer, autumn, and winter) to determine its Weibull parameters and the characteristics of wind speed. To achieve this, a statistical analysis was conducted, and with this information, wind can be utilized effectively as a resource for electric generation.

The main objective of this study was to identify seasonal wind characteristics to assess them and determine their potential using different types of wind turbines based on their power curves and power coefficients. We considered a very important step in the process of wind turbine selection, which is the assessment or characterization of the wind speed. In this study, we also include a proposal to relate the wind turbine efficiency through its power coefficient and the conditions of the wind at each specific site.

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

#### *2.1. Sites and Data*

Mexico has five states along the Gulf of Mexico: Tamaulipas, Veracruz, Tabasco, Campeche, and Yucatan, as seen in Figure 1.

In order to carry out this study, data for 141 different sites along the Mexican Gulf were collected from the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) from the National Aeronautics and Space Administration (NASA); these are re-analysis data that are long-term, model-based analyses of multiple datasets using a fixed assimilation system [32]. MERRA-2 has presented data every hour from 1980 until 2017. According to Shang [33] the network measures derived from the empirical observations are often poor estimators of the true structure of system as it is impossible to observe all components and all interactions in many real-world complex systems. This problem occurs when there is missing data; in this study, MERRA-2 has no missing data; therefore, this problem is avoided and can be considered to be unbiased data.

In Figure 2, we give the geographic positions where MERRA-2 data were taken from.

**Figure 1.** Mexican states along the Gulf of Mexico.

**Figure 2.** Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) data locations.

#### *2.2. The Weibull Distribution and Seasonal Wind*

There are some statistical distribution functions for analyzing wind data, such as the lognormal, normal, Rayleigh, and Weibull probability distributions [34,35]. The Weibull function is the most-used function to assess wind energy potential because shows variables as shape and scale parameters [36], where these parameters are obtained using estimation methods, such as the maximum likelihood method, and the goodness of the resulting fits are evaluated using several indicators, e.g., the coefficient of determination (R2) [24,37–39]. The R2 is a statistical measure that represents the proportion of the variance for a dependent variable that is explained by an independent variable; R<sup>2</sup> is generally interpreted as the percentage of a value's movements that can be explained by movements of another variable [40].

#### 2.2.1. Wind Model

The Weibull distribution and cumulative distribution functions are expressed in Equations (1) and (2), respectively [41]:

$$f(v) \;= \frac{k}{c} \Big(\frac{v}{c}\Big)^{k-1} \exp\left[-\left(\frac{v}{c}\right)^k\right] (k > 0, v > 0, c > 1),\tag{1}$$

$$F(v) = 1 - \exp\left[\left(\frac{v}{c}\right)^k\right] \tag{2}$$

where *k* is the shape parameter (dimensionless), which is considered a Weibull form parameter because it specifies the shape of the distribution taking place within values between 1 and 3. A small value for *k* signifies very variable winds, while constant winds are characterized by a larger *k* [42]; when the shape parameter is 2, it is considered to represent Rayleigh distribution. The scale parameter *c* has the same units as wind speed (m/s) and is proportional to the mean wind speed (*vm*), where *v* (m/s) is the wind speed registered in the site.

If the mean wind speed in Equation (3), and the standard deviation (σ) are known, *k* and *c* can be calculated using Equations (4) and (5) as follows [43]:

$$\left. \psi\_m \right| = \left. \frac{1}{n} \right| \sum\_{i=-1}^{n} \psi\_i \Big| \tag{3}$$

$$k = \left(\frac{\sigma}{\upsilon\_m}\right)^{-1.086} \text{ .} \tag{4}$$

$$x = \frac{v\_m}{\Gamma\left(1 + \frac{1}{k}\right)}.\tag{5}$$

The gamma function, Γ, can be calculated using Equation (6):

$$
\Gamma(r) = \int\_0^\infty x^{r-1} e^{-x} dx. \tag{6}
$$

#### 2.2.2. Wind Power Analysis

The observed wind power density (*WPDO*) can be obtained using Equation (7):

$$WPD\_O = \frac{\sum\_{i=1}^{n} 0.5 \rho v\_i^3}{n} \,\tag{7}$$

where ρ is the air density and is calculated using the ideal gas law (see Equation (8)), in which *T* is the absolute temperature (K), *p* is the absolute pressure (Pa), and *R* is the specific gas constant (J/kg·K): As an approximation, *R* = 0.286 is used for dry air.

$$
\rho = \frac{p}{RT} \tag{8}
$$

#### *2.3. Wind Variables Analysis*

Calculations for wind variables were performed by implementing a self-written code on the R free software environment and language [44]. This environment, besides being free, provides powerful tools for statistical computing and graphical display via eight packages, and when required, they can be extended with additional packages available through the comprehensive R archive network (CRAN) family available on the Internet. In fact, in the former study, six base packages were used (tools, stats, graphics, grDevices, utils, and base) [44], which were complemented with the extra packages (lubridate [45], RColorBrewer [46], ggplot2 [47], and gridExtra [48]).

#### 2.3.1. Data Processing

An iterative process was run for all 141 studied sites, which were distributed in a grid-shaped arrangement along the coast of the Gulf of Mexico (see Figure 2) using the previously mentioned MERRA-2 files as input data. These files contained complete time series data (no missing values) of the surface pressure (kPa), air temperature (◦C) at 2 m and 10 m above ground level, wind speed (m/s) at 50 m above ground level, and wind direction (◦) at 60 m above ground level; all these variables have hourly records from 1980 to 2016.

An annual and a seasonal analysis was performed for every site, where for the sake of simplicity, spring was considered to include the months of March, April, and May; summer included June, July, and August; autumn included September, October, and November; and Winter included December, January, and February.

The structure of the program consisted of an initial pre-processing of the MERRA-2 data, which generated an output file that was used for obtaining the geolocated values of the wind variables, as well as customized plots via two scripts named Weibull Analysis and Directional Analysis, as can be observed on Figure 3, and is explained as follows.

**Figure 3.** Program structure. CDF: cumulative density function, PDF: probability distribution function.

#### 2.3.2. Pre-Processing

In the pre-processing step, a time adjustment was performed on the MERRA-2 data in order to assign the corresponding UTM zone (for Mexico −6 h) and make summertime corrections. Data was classified for years, months, days, and hours for further manipulation, and the air density and observed WPD were obtained via Equations (9) and (10), respectively.

Once the pre-processing was finished, an output was generated for feeding the following two independent and complementary scripts.

#### 2.3.3. Weibull Analysis

The Weibull probability distribution function (PDF) and cumulative distribution function (CDF) were obtained and plotted according to Equations (1) and (2), respectively, for the annual and seasonal series.

Tables were generated with Weibull shapes and scale factors (*k*, *c*), mean wind speed (*vm*), its corresponding standard deviation (σ), Weibull most-probable wind speed (*vmp*), mean air density (ρ*m*), and the observed wind power density (*WPDO*).

#### 2.3.4. Directional Analysis

In order to complement the Weibull distribution, a directional analysis was performed. For this, data was grouped according to the wind direction into 12 bins, corresponding to sectors as 30◦ segments. The cumulative *WPDO* was obtained for every sector by considering the one with the highest value as the prevailing direction, which was obtained for complete years, as well as for every season.

#### *2.4. Wind Clustering*

The wind performance was grouped into several clusters using a *k*-means clustering model to identify its monthly behavior, which was analyzed to diagnose the time where the wind speed could be used to generated power. The *k*-means algorithm is widely used in data mining for the partitioning of *n* measured quantities into *k* clusters [49]; according to Sugar and James [50], the classification of observations into groups requires computing the distance between the observations [51–54]. The *k*-means algorithm is one of the simplest unsupervised machine learning algorithms, where unsupervised algorithms make inferences from datasets using only an input vector without referring to known, or labelled, outcomes [53]. We can define a cluster as a data set with similar characteristics. The *k*-means algorithm identifies *k* centroids and then allocates every data point to the nearest cluster [53]. This algorithm has been used in a study done by Wang et al. [55], where the *k*-means clustering algorithm was used to find the largest historical samples that had the greatest influence on forecasting accuracy to improve the efficiency of the proposed model.

A method for clustering was proposed by Deng et al. [52], which used a Weibull distribution to establish that an unclustered dataset *P* can be represented using Equation (9):

$$P = \begin{cases} p\_{i\prime} \text{ i } = \text{ 1, 2, 3, \dots, N}^{i} \text{ } \end{cases} \tag{9}$$

where *pi* denotes the *i*th element and *N<sup>i</sup>* is the number of observations. The set of clusters *C* can be given as:

$$\mathcal{C} = \left\{ \mathbb{C}\_{\bar{j}} | \mathbb{C}\_{\bar{j}} \subset P, \ j = 1, 2, 3, \dots, N^{\bar{j}} \right\} \tag{10}$$

and the set of observations within a cluster, *Cj*, is represented by:

$$\mathcal{C}\_{\vec{\lambda}} = \left\{ \mathcal{C}\_{\vec{\lambda}} \middle| \mathcal{C}\_{\vec{\lambda}} \subset P, \ j = 1, 2, 3, \dots, N^{\vec{\lambda}} \right\}. \tag{11}$$

*Cjk* and *Njk* are the *k*th observation with the *j*th cluster and the number of observations, respectively, in cluster set *Cj*. The set of centroids associated with the clusters *W* is given as:

$$\mathcal{W} = \left\{ \mathcal{W}\_{\dot{j}} \middle| \mathcal{W}\_{\dot{j}} \subset \mathcal{C}\_{\dot{j}}, \ j = 1, 2, 3, \dots, N^{\dot{j}} \right\}. \tag{12}$$

where *Wj* is the *j*th centroid.

For wind clustering, each observation *P* is assigned to a cluster *Cj* and its centroid will be represented by its mean wind speed.

#### *2.5. Wind Turbine Selection*

The energy available for conversion mainly depends on the wind speed and the swept area of the wind turbine. Using Newton's Law *F* = *ma*:

$$E = \text{mas}\_{\text{\textquotedblleft}}\tag{13}$$

where *E* is the kinetic energy, *m* is the mass, *a* is the constant acceleration, and *s* is the distance. From kinematics, the acceleration can be expressed using Equation (14):

$$a = \frac{\left(v^2 - u^2\right)}{2s}.\tag{14}$$

Because the initial velocity of the object is 0, then *u* = 0, such that:

$$a = \frac{v^2}{2s}.\tag{15}$$

Substituting Equation (15) in Equation (13) gives:

$$E = \frac{1}{2}mv^2.\tag{16}$$

The rate of change of energy of the power in the wind is given by:

$$P = \frac{dE}{dt} = \frac{1}{2}v^2 \frac{dm}{dt}.\tag{17}$$

The mass flow rate is expressed as:

$$\frac{dm}{dt} = \rho A \frac{dx}{dt'}\tag{18}$$

and the distance's rate of change is:

$$\frac{d\mathbf{x}}{dt} = v.\tag{19}$$

Substituting Equation (18) into Equation (19) gives:

$$\frac{dm}{dt} = \rho Av.\tag{20}$$

Hence, from Equation (17), the wind power output generated by a wind turbine can be expressed using Equation (21):

$$P\_{WT} = \frac{1}{2} \rho A l I^3,\tag{21}$$

where *PWT* is the rated power, ρ is the air density, *A* is the rotor area, and *U* is the wind speed approaching the turbine. According to Grillo et al. [56], the power coefficient is a function of the tip speed ratio, known as lambda (λ), and the blade pitch angle (β) (◦). λ can be calculated Equation (22):

$$
\lambda = \frac{a\mathbb{R}}{\mathcal{U}}\!\!/\tag{22}
$$

where ω is the rotational speed of the rotor (rad/s) and *R* is the rotor radius (m).

To characterize the wind speed in this study, a proposal for wind turbine selection is given, where this proposal considers the air density, rotor area, wind turbine power curve, and power coefficient (*Cp*).

Song et al. [57] defined the power coefficient *Cp* as the ratio of the power extracted from the wind turbine to the available power. *Cp* can be calculated as follows:

$$\mathbb{C}\_p = \frac{P\_n}{\frac{1}{2}\rho A u\_j^3}.\tag{23}$$

The *Cp* of a wind turbine is a measurement of how efficiently the wind turbine converts the energy in the wind into electricity.

Wind speed is one of the most important parameters in determining the electric power, and the general equation is related to the density of air, wind speed, and swept area, as in Equation (21), and represents the total energy obtained from the wind resource; however, in terms of generating electricity, only a certain proportion of energy can be converted and is expressed by Equation (24), as follows:

$$P\_{\mathcal{E}} = \eta\_{\mathcal{E}} \eta\_{\mathcal{W}} C\_p P\_{\mathcal{W}\mathcal{T}\prime} \tag{24}$$

where *Pe* is the amount of electric power generated, η*<sup>e</sup>* is the electrical conversion efficiency of the wind turbine, and η*<sup>m</sup>* is the mechanical efficiency [58].

There is an optimization proposal that uses the wind turbine efficiency as a fundamental variable to determine the power output generated by a wind power farm [59]. In this proposal, as in this study to determine it, a Weibull distribution was used. The parameter *c* from the Weibull distribution was used to select a wind turbine according to its *Cp*; in this case, the maximum efficiency of a wind turbine was compared to the parameter *c* calculated from the wind speed site studied.

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

The annual Weibull parameters *c* and *k*, *WPDF*, *WPDO*, ρ, *vm*, *vmp*, and *vmaxE* were obtained for the 141 MERRA-2 data locations. Table 1 only presents the places with the highest and lowest mean wind speeds for each state (colored blue in Figure 2); these values were used as references for the seasonal analysis.


**Table 1.** MERRA-2 analysis.

Tamaulipas was the state with more zones with high values of wind than the other states, where its maximum *vm* = 7.34 m/s, scale and shape parameters were between 4.17–8.26 m/s and 1.97–2.91, respectively. In this case, by comparing *vm*, *c*, and *k*, it can be established that in the time domain, most of the wind resource had values higher than its mean wind speed. In Veracruz, the scale and shape parameters showed values between 3.86–6.88 m/s and 1.72–2.41, respectively. Tabasco presented more wind speed and wind resource than Veracruz with Weibull parameters ranging between 4.02–7.48 m/s for *c*, and 2.51–3.07 for *k*. Campeche and Yucatan showed a higher *k*'s than the others states, with 2.71 and 2.95, respectively, which can be interpreted as a long period of time with winds higher than their means; *c* was 7.47 m/s and 7.88 m/s for Campeche and Yucatan, respectively.

According to Katinas et al. [29], the Weibull parameter *c* has the same behavior as the WPD, where in all cases, when the scale factor grows, the WPD grows as well. It is important to mention that the lowest shape parameter corresponds to the lowest wind speed, which means that the frequency of data below the mean was greater than that above it.

#### *3.1. Seasonal Weibull Parameters and WPD*

Seasonal wind speed characterization was carried out by dividing the data according to the four annual seasons: spring (March, April, and May); summer (June, July, and August); autumn (September, October and November); and winter (December, January and February). For the sake of analysis, two sites for each state were chosen by considering its highest and lowest mean wind speed (see Figures 4–8).

In Figure 4, it can be observed that the points Tam34 and Tam48 had the highest and the lowest mean wind speed—8.26 m/s and 4.18 m/s, respectively—in the state of Tamaulipas. At the annual level, Tam34 had *k* = 2.71, which represents a very good value due its higher magnitudes of wind speed; at the seasonal level, the Weibull distribution shows that spring and winter are the periods of time during an average year where there was more available wind resource.

In the state of Veracruz, 31 locations were studied, with the highest scale parameter being found for Ver29, where this was located is in the Isthmus of Tehuantepec, one of the windiest zones in the world. During an average year, its highest mean wind speed was 5.72 m/s, its *c* was 6.46 m/s, and *k* was 2.40, and had a windy season during spring, autumn, and winter (see Figure 5).

In Tabasco (Figure 6), Tab2 had *c* = 6.90 m/s, *k* = 2.84, and a mean wind speed of 6.15 m/s; these values represent a good wind resource because analyzing Weibull parameters showed that *c* was higher than its mean, and *k* showed that most of the time, the wind speed data was above its mean. In contrast, the seasonal Weibull of Tab9, which was the point with the lowest values in this state (*c* = 4.02 m/s, *k* = 2.76, *vm* = 3.58 m/s), allowed to determine a period of wind during an average year that had constant magnitudes of wind values, namely winter, spring, and summer, as shown.

The analysis of Campeche (Figure 7) showed its windiest location was Cam14, which had a mean wind speed of 5.82 m/s, and *c* and *k* values equal to 6.51 m/s and 3.08, respectively. The lowest valued location, Cam20, had *c* = 2.90 m/s and *vm* = 2.64 m/s, where both values are similar with a difference of 0.26 m/s, which can be considered as a location that did not have variations during an average year. Its *k* = 3.04 could be considered a very good frequency of wind speed; however, due its low wind speed, it did not represent an impactful wind resource, which can be seen in the seasonal analysis because its variation was minimal.

In Yucatan, as seen in Figure 8, its seasonal wind variation was divided in two periods, highest one between winter and spring and the lowest one between summer and autumn. The location with the highest wind speed was Yuc2 (*vm* = 4.56 m/s, *c* = 5.10 m/s, and *k* = 3.07), and the lowest location, Yuc17, was inside a rainforest, which reduced the magnitude of the wind speed (*vm* = 2.67 m/s, *c* = 2.99, and *k* = 3.10).

**Figure 4.** Annual and seasonal Weibull distributions in the state of Tamaulipas.

**Figure 5.** Annual and seasonal Weibull distributions for the state of Veracruz.

**Figure 6.** Annual and seasonal Weibull distributions for the state of Tabasco.

**Figure 7.** Annual and seasonal Weibull distributions for the state of Campeche.

**Figure 8.** Annual and seasonal Weibull distributions for the state of Yucatan.

These results are in agreement with the study developed by Herrero-Novoa [30], where it was established that dividing an average year into two periods (spring–summer and fall–winter) provides more accuracy regarding wind characterization; in a similar fashion, for the current study, an average year was divided into four periods of time (seasons) where even more accurate wind periods were established. This allowed us to determine the wind characteristics, specifically its WPD, in the most relevant period. By observing Figure 8 and considering the values of the parameters obtained in Table 1, we established that in Yuc2, the parameter *c* = 5.10 m/s was higher than its average of 4.56 m/s, and its *k* = 3.07; therefore, at this site, the wind speed was over 4.56 m/s most of the time. Furthermore, according to Vazquez et al. [42], this value of *k* means that the site had constant winds.

In addition to the Weibull analysis, in Figures 9 and 10, the prevailing wind directions and the magnitude of its *WPDO* can be observed for the annual and seasonal behavior, respectively.

**Figure 9.** Annual prevailing wind direction and its *WPDO* along the states of the Gulf of Mexico.

The highest *WPDO* points were along Tamaulipas and Veracruz, although Tamaulipas had more locations than Veracruz, while Veracruz has the highest resource next to the sea, Tamaulipas had it at the border with the United States of America.

**Figure 10.** Seasonal prevailing wind direction and its *WPDO* along the states of the Gulf of Mexico.

As can be observed in Figure 10, the seasons with greatest *WPDO* were autumn and winter, and the states with greatest resource were Tamaulipas and Veracruz. Regarding wind orientation, this figure shows that Tamaulipas had two periods with different prevailing directions, with the first one during spring and summer and the second one in autumn and winter. The WPD in southern Veracruz presented the same direction most of the time (southeast), meanwhile Tabasco, which is the southernmost state along the Gulf of Mexico, had different prevailing orientations throughout an average year. It can also be seen that Campeche had a prevailing direction in winter and spring, and for Yucatan, autumn, winter, and spring had northeast as the predominant direction, while it was southeast in summer.

#### *3.2. Wind Speed and Wind Power Clustering*

Using *k*-means clustering, the wind speed could be clustered into months, and in these clusters, the power output could be calculated. Figure 11 shows the clustering done for wind speed from 1980–2017.

**Figure 11.** Wind speed clustering in the Mexican states of the Gulf of Mexico.

Wind speed was divided into three clusters—C1, C2, and C3—as shown in Figure 11 where the blocks represent the amount of wind speed each month.

Taking these results and knowing the windy seasons throughout an average year, the wind power output could be calculated. A critical point in this stage was to select the wind turbine, where in Table 2, we present the *Cp* of different wind turbines as a first step to know the efficiency of each one of them. The wind turbines were selected according to parameter *c* of the points analyzed.



Six wind turbines were found such that their *Cp* fit with the parameter *c* calculated from the wind speed data. Table 3 presents the wind turbine selected for each point along the Gulf of Mexico.


**Table 3.** Wind turbine selection.


**Table 3.** *Cont.*

A total of 58 points were fitted with an appropriate wind turbine *Cp*, where the wind speed of the site was related to the wind speed at the maximum wind turbine efficiency; with this information, the best option for the wind turbine to be used could be determined based on its *Cp*. The probable energy production was calculated for each wind turbine selected using Equation (16) based on the electrical and mechanical efficiency given by the manufacturer.

#### **4. Conclusions**

The seasonal characterization of wind speeds along the Mexican states of the Gulf of Mexico was done for 141 locations with MERRA-2 data, with records between 1980–2017. An average year of wind data for each location was obtained and these were divided according to the seasons of the year (spring, summer, autumn, and winter). This distinction allowed us to describe the wind characteristics more accurately, especially the WPD, which allowed for establishing wind seasons properly, as well as the description of the prevailing wind direction. It is was found that there were different wind seasons along the states of the Gulf of Mexico and they were established for each state. The wind season in Tamaulipas occurred in winter and spring; Veracruz and Tabasco had a wind season in autumn, winter, and spring; and in Campeche and Yucatan, the wind season was in winter and spring, where these states shared their coast with both the Gulf of Mexico and the Caribbean Sea. This variation allowed us to determine that the wind season could depend upon other factors regarding the climatology of each state.

The state with the greatest wind resource was Tamaulipas, as seen in Figure 9; it had plenty of locations with good values for wind, scale, and shape parameters, where its highest *vm*, *c*, and *k* was 7.34 m/s, 8.26 m/s, and 2.91 respectively. North and south Veracruz had the highest values of wind speed, where its southern zone corresponded to the Isthmus of Tehuantepec, one of the zones with the greatest wind energy resource in the world; its highest parameters were *vm* = 6.10 m/s, *c* = 6.88 m/s, and *k* = 2.41. Tabasco presented the following parameters: *vm* = 6.67 m/s, *c* = 7.48 m/s, and *k* = 3.07. These three states (Tamaulipas, Veracruz, and Tabasco) presented the same wind season.

The highest parameters for Campeche were *vm* = 6.57 m/s, *c* = 7.39 m/s, and *k* = 2.71, and the Yucatan parameters were *vm* = 7.88 m/s, *c* = 7.04 m/s, and *k* = 3.07; these two states presented the same wind season.

Three clusters were found, which were useful for determining the wind speed in monthly groups and to visualize the wind stations in the study area of the Gulf of Mexico.

Using *Cp*, optimal wind turbines were determined by comparing it with the parameter *c* of Weibull distribution. In this study, 58 sites were fitted with a wind turbine *Cp*, and by analyzing these data, three wind turbines were selected: Acciona AW70/1500 Class III with a *Cp* = 0.4188 at 6.5 m/s, Acciona AW70/1500 Class II with a *Cp* = 0.4457 at 7.5 m/s, and Vestas V90 1.8 MW with a *Cp* = 0.4381 at 8 m/s.

At the end of this study, the probable energy production was calculated for each wind turbine applied at each site with its own respective conditions.

Regarding future research lines, we propose following the proposal to select the most efficient wind turbine by adding more constraints, such as wind direction, roughness, and orography.

**Author Contributions:** Conceptualization, A.-J.P.-M., G.A. and Q.H.-E.; methodology, A.-J.P.-M., G.A. and Q.H.-E.; formal analysis, A.-J.P.-M., G.A. and Q.H.-E.; investigation, A.-J.P.-M., G.A. and Q.H.-E.; resources, A.-J.P.-M., G.A. and Q.H.-E.; writing—original draft preparation, A.-J.P.-M., G.A. and Q.H.-E. All authors have read and agreed to the published version of the manuscript.

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

**Acknowledgments:** The authors acknowledge the Mexican "Secretaria de Educacion Publica" for the support to the research group UV-CA-466 "Mecanica Electrica" in the project "Evaluación del potencial eólico para disminuir el pico de demanda eléctrica en zonas tropicales. Caso Coatzacoalcos Veracruz México.".

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

#### **References**


© 2019 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* **Energy Recovery from Waste Tires Using Pyrolysis: Palestine as Case of Study**

**Ramez Abdallah 1, Adel Juaidi 1, Mahmoud Assad 1, Tareq Salameh <sup>2</sup> and Francisco Manzano-Agugliaro 3,\***


Received: 22 February 2020; Accepted: 3 April 2020; Published: 9 April 2020

**Abstract:** The first industrial-scale pyrolysis plant for solid tire wastes has been installed in Jenin, northern of the West Bank in Palestine, to dispose of the enormous solid tire wastes in the north of West Bank. The disposable process is an environmentally friendly process and it converts tires into useful products, which could reduce the fuel crisis in Palestine. The gravimetric analysis of tire waste pyrolysis products from the pyrolysis plant working at the optimum conditions is: tire pyrolysis oil (TPO): 45%, pyrolysis carbon black (PCB): 35%, pyrolysis gas (Pyro-Gas): 10% and steel wire: 10%. These results are depending on the tire type and size. It has been found that the produced pyrolysis oil has a High Heating Value (HHV), with a range of 42 − 43 (MJ/kg), which could make it useful as a replacement for conventional liquid fuels. The main disadvantage of using the TPO as fuel is its strong acrid smell and its low flash point, as compared with the other conventional liquid fuels. The produced pyrolysis carbon black also has a High Heating Value (HHV) of about 29 (MJ/kg), which could also encourage its usage as a solid fuel. Carbon black could also be used as activated carbon, printers' ink, etc. The pyrolysis gas (Pyro-Gas) obtained from waste tires mainly consist of light hydrocarbons. The concentration of H2 has a range of 30% to 40% in volume and it has a high calorific value (approximately 31 MJ/m3), which can meet the process requirement of energy. On the other hand, it is necessary to clean gas before the burning process to remove H2S from Pyro-Gas, and hence, reduce the acid rain problem. However, for the current plant, some recommendations should be followed for more comfortable operation and safer environment work conditions.

**Keywords:** pyrolysis; solid tire wastes; PCB; TPO; Pyro-Gas; industrial scale; Palestine

#### **1. Introduction**

The world is facing a huge energy gap due to the increase of the human population and the unexpected growth in industry around the world. In Palestine, the consumption of energy has hugely increased in the last years. According to the Palestinian Central Bureau of Statistics, Palestinians consumed more than 30 to 40 thousand barrels per day of petroleum in 2017 [1]. However, Palestine suffers from a scarcity of natural resources and mineral wealth. The suffering of the Palestinians mostly lies in the scarcity of conventional energy sources, such as oil and gas, and additionally, the expensive price of energy, which is equivalent to the most expensive cities in the world. Inspired by this suffering, Palestinians seek to overcome this dilemma and to find alternative sources for conventional fuels.

Nowadays, the trend in the technology of energy is heading towards alternative sources and mainly, renewable energy sources, in order to overcome the rising demands of energy, which are essential in diverse fields of life. There are many sources for alternative fuels, which are based on recycling industrial and household wastes. These wastes exist in different forms: liquid wastes (e.g., cooking oil) and solid wastes. The solid wastes are any kind of rubbish, tires, paper, cardboard, plastic, wood, glass, metallic wastes and ashes.

The world in general, and Palestine in particular, are facing the problem of solid waste due to the highly increasing numbers of the population, and the lack of material resources needed to solve this problem and reduce its risks. Therefore, the solid waste is randomly collected and assembled in the dumps, contrary to the conditions of environmental health, but this is still not a suitable solution. The waste to energy process is a very simple concept involving recovering the maximum energy from waste with minimum damage to the environment by different thermal processes, such as the pyrolysis process [2]. With the rapid development of the auto industry, "black pollution" is being created by more and more waste tires, which now threaten the living environment of humans.

Millions of tons of waste tires are produced worldwide: about 2.5 million tons is generated in North America, Japan generates one million tons and the European Community generates more than 2.5 million tons [3]. In Palestine, in 2012, it was estimated that more than 380 thousand scrap tires were generated from different sources [4]. The tires are composed of vulcanized rubbers, which include natural rubber (NR), styrene butadiene (SBR), and polybutadiene (BR), carbon black, steel wire, cords and additives.

Getting rid of non-biodegradable solid waste tires is one of the most important environmental problems in the world. The durability and immunity of these wastes towards biological degradation makes disposal and re-treatment difficult [2,5]. In Palestine, the solid waste tires are randomly dumped in landfills. These landfills usually face accidental fires that produce extremely thick black smoke and toxic emissions. Such emissions are deleterious to health and may cause the acid rain [5]. Therefore, the biomass pyrolysis is a suitable solution to recycle these waste tires with more environmentally friendly methods. Pyrolysis is considered a cheap way to obtain alternative liquid fuel. By 2050, the International Energy Agency expects that 27% of the fuel used in the transportation industry may be sustainably produced from biomass and waste resources [6].

Tire recycling is the conversion of scrap rubber tires into hydrocarbons for reuse as oils and fuels, and the recovery of steel and carbon is one of the more desirable methods for economically reclaiming scrap tires. Tire recycling is achieved by the pyrolysis process, see Figure 1.

Pyrolysis is the thermochemical decomposition of tires at high temperatures without oxygen. Pyrolysis takes place under low pressure and high temperatures above 420 ◦C. Practically, pyrolysis cannot be performed oxygen free, since with a limited quantity of oxygen present in the reactor, oxidation will occur at a limited quantity. Pyrolysis of waste tires will produce pyrolysis black carbon (PCB) (30–35 wt%), crude tire pyrolysis oil (TPO) (40–45 wt%), scrap steel (10–15 wt%) and pyrolysis gas (Pyro-Gas) (10–15 wt%) [7].

**Figure 1.** Products of the tire pyrolysis process.

#### *Pyrolysis Technology and Energy Security in Palestine*

Lately, increasing energy consumption, high-price volatility in fossil fuels, interruptions in energy imports and exports, climate change and pollution of the environment are just a few of the many energy security issues in countries throughout the world [8]. Renewable energy and bioenergy play a substantial role in economy, social and environmental life [9,10], bringing significant benefits for developing countries, specifically Palestine. Besides, use of these energy sources avoids or eliminates the negative environmental impacts. Regarding the political situation in Palestine and to avoid the highly frequent blackouts, use of renewable energy and bioenergy (such as pyrolysis) will be a sustainable solution to guarantee the energy supply in the long term. A significant driver of pyrolysis is also the growing interest in environmental sustainability [11]. Energy security is a significant factor that enables the system to function optimally and sustainably, free of threats and risks [12]. Energy security guarantees energy supply reliability at an affordable price, and is therefore affected by economic, political and technical conditions [13].

Energy supply security and environmental issues are considered to be among the main concerns for most countries [8]. While Palestine is highly dependent on imported energy, the results show that solar energy and bioenergy (Pyrolysis) are the most promising renewable and sustainable energy sources. The advantages of pyrolysis include a highly flexible technology that can handle a wide variety of feedstocks consisting of forest residues, food waste, tire waste and solid municipal waste [11]. Palestine is among the most promising locations for utilization of bioenergy and renewable energy sources because of the geographic location of this area, which supplies it with plenty of high solar radiation, and the availability of bioenergy sources. Also, it is highly dependent on energy imported from Israel [9].

To achieve sustainability, this project involves producing a renewable fuel from the pyrolysis process of waste tires. The fuel production from this process involves less costs and less greenhouse gas emissions and makes the environment cleaner. The sustainability aims of this project are directed at solving three main issues: economic, environment and social, in order to achieve energy security and a clean environment in Palestine.

The present study aims at solving the firm technical problems of tire wastes in Palestine. The novelty of this research is that it adds sustainability as an aspect of energy security in Palestine's energy systems. This is done by introducing the first commercial plant of tire pyrolysis in Palestine.

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

The disposal of waste tires is considered a main economic and environmental problem in the waste management field. For the time being, attention is growing in pyrolysis due to its ability to make products of economic value in different fractions by changing the operating conditions, such as the operating temperatures and heating rate [14]. Pyrolysis also transforms waste tires into tire pyrolysis oil (TPO), pyrolysis carbon black (PCB) and pyrolysis gas (Pyro-Gas), which are considered as high value fuels and chemicals [15,16]. Table 1 below summarizes the pyrolysis research during the period from 1980 to 2018.



#### *Energies* **2020**, *13*, 1817


Pyrolysis Gas: Pyro-Gas. The pyrolysis process allows the communities to have an alternative environmentally friendly source of energy, reduce the need for non-renewable fossil minimize the landfilling problems.

and

#### *Energies* **2020**, *13*, 1817

#### **3. Methodology**

As mentioned earlier, this study was conducted in Palestine to solve the environmental problem of direct burning of the tires and reduce pollution. It is believed that the pyrolysis process of such waste is a suitable solution. The classical tires consist of different materials such as synthetic rubber, steel wires and carbon black. These materials are homogenously distributed over the tire volume. To dispose of these materials using pyrolysis technology, the tires are usually shredded to smaller pieces, washed with water and dried. The tires are placed in a rotary reactor to extract the wires and produce gaseous and liquefied oil. Since no oxygen exists inside the reactor, the carbon black and steel properties will not change.

This study investigates the pyrolysis of tires using an industrial-scale system. The industrial scale is basically a commercial pyrolysis plant, on about 6000 m2 of land in the north of Jenin City, West Bank, Palestine. The pyrolysis plant consists of five sections: the main plant section, the offices section, the laborers rooms section, the waste tire section and the warehouse for the products from pyrolysis. This plant was established in 2013 to get rid of the enormous tire wastes in the north of West Bank. The flow diagram of the pyrolysis process system in this plant is shown in Figure 2.

**Figure 2.** Flow diagram for the tire pyrolysis process.

The plant consists of a horizontal axis cylindrical reactor (6 m long and 2 m diameter) that has a batch mode and rotates with a rotational speed of 0.4 rpm. The reactor chamber can fit up to five tons of partly shredded tires (only small side rings are removed, and steel wires are kept inside). The cylindrical chamber is made from steel and it is initially heated by flames from external source power. Once the system starts producing gaseous fuel from the pyrolysis process, the chamber is now heated by burning the produced gaseous fuel. Therefore, the system is self-sustained in terms of energy. Electricity was used to rotate the chamber for better stirring. The tires were added, and a heavy door was closed tightly. The heat exchangers consist of two groups. The first group contains three heat exchangers that condensate the passing gas coming from the separator. The second group contains two heat exchangers to condensate the remaining gas coming out from the first group. The condensate fuel is filtered to purify the fuel from the contaminations (such as carbon black).

The heating system is composed of an insulated furnace below the reactor. The furnace contains two liquid burners and two gas burners. An insulating cover surrounds the reactor drum. The resulting hot gases pass between the reactor drum and its insulated cover and give enough heat for the pyrolysis process.

Two storage tanks were used in this system: the first tank stores the produced fuel from the first group of heat exchangers, while the second tank stores the remaining produced liquid fuel. The fuel in the first tank is heavier (higher molecular weight fractions) and less clear than the fuel in the second tank. However, the heavy fuel is still demanded in the market for its high heating value and efficient burning. The excess gases that do not liquefy are stored into a third tank, which feeds the heating process. The exhaust gas from the combustion process is filtered for toxic gases before it is spilled out to the atmosphere. The reactor works at a relatively low pressure (the maximum reactor pressure is 0.08 bars). Figure 3 shows the reactor temperatures during heating and cooling processes for four different runs. The total tire load for each run (batch) was 5000 kg. Firstly, the pyrolysis oil produced from previous pyrolysis runs is burned to heat the reactor for the first 3 to 4 h until the temperature reaches 420–500 ◦C. Secondly, the temperature is held at 420–500 ◦C by burning the Pyro-Gas produced from the current pyrolysis run. Finally, the reactor is left to cool down naturally for around 10 h. When the reactor is cooled below 60 ◦C, steel wire and pyrolysis black carbon (PCB) are removed. The PCB fraction is in the form of powder and separated from steel wire while the reactor drum is continuously rotating.

**Figure 3.** The reactor temperatures during heating and cooling for the four experiments.

The flue exhaust gas system consists of a chimney, a desulfurization and dust removal unit and an exhaust suction fan. The flue gas is passed through the desulfurization and dust removal unit to filter the dust and remove the sulfur from the flue gas. Finally, the filtered flue gas is exhausted to the atmosphere.

#### *Analaytical Techniques*

The original sample of tire and the sample of pyrolysis carbon black (PCB) for the four experiments underwent thermogravimetric analysis to take simultaneous results of differential thermogravimetric analysis (DTA) signals. The thermogravimetric analyses are done by burning the samples (approximately 15 mg for each sample) with a constant heating rate and gas flow rate (15 ◦C/min) and (100 mL/min), respectively. The calorific value, proximate and elemental analyses were performed, using the procedures of American Society for Testing and Materials (ASTM) standards.

The four samples of tire pyrolysis oil (TPO) from different experiments in the industrial pyrolysis plant were homogenized by mixing the oil well before the analysis was carried out. ASTM standard procedures were also used to find the physical properties of tire pyrolysis oil (TPO): flash point, viscosity, pour point, density and HHV. Elemental analysis of TPO was performed by an elementary analyzer.

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

#### *4.1. Yields of the Pyrolysis Products*

The yields for each of the fractions of the pyrolysis process obtained from the four experiments are shown in Table 2. The total tire load for each experiment (batch) was 5000 kg. The products of the pyrolysis process for the four experiments were observed to be similar with small differences. These differences are attributed to the tires' diversity and their size.


**Table 2.** Yield of pyrolysis products.

#### *4.2. Characterization of Solid Tires and Pyrolysis Carbon Black (PCB)*

The original tire and the carbon black fraction for the four experiments were subjected to thermogravimetric analyses. Calorific value and elemental analyses were done according to the ASTM standard procedure.

Proximate and elemental analyses were done on the original tire and the pyrolysis carbon black (PCB) fractions obtained from the four experiments and the results are tabulated in Table 3. It can be concluded that small differences were obtained between the PCB fractions for the four experiments. It can also be noticed that chlorine and sulfur were condensed in the PCB fraction. Approximately, most of the chlorine existing in the waste tires was concentrated in the PCB fraction. Therefore, it is important to take this into consideration once PCB is used as a fuel. So, a system for cleaning the gas would be necessary for the residual gases resulting from combustion of PCB. The PCB residue of the pyrolysis tire process will produce three times the amount of chlorine produced by the original tires [40]. The High Heating Value (HHV) for pyrolysis carbon black (PCB) is about 36 MJ/kg. On the other hand, the HHV for convectional liquid petroleum is 39.5 to 45 MJ/kg. The HHV of PCB is about two times that of wooden biomass (18.5 to 21 MJ/kg). From the analysis results above, it could be concluded that PCB is a possible source of energy.

Figure 4 shows a Derivative Thermogravimetry (DTG) for the original tire and the pyrolysis carbon black (PCB) fractions obtained from the four experiments. It can be concluded from Figure 4 that the highest mass loss rate for the original tire occurs nearly at 390 ◦C, while for the PCB, it occurs at 570 ◦C.


**Table 3.** Proximate, elemental and heating value analysis for the tire and carbon black.

<sup>1</sup> Results in dry basis and ash free.

**Figure 4.** Derivative Thermogravimetry (DTG) plots of the original tire and the pyrolysis carbon black (PCB) fractions obtained from the four experiments.

#### *4.3. Characterization of Tire Pyrolysis Oil (TPO)*

The tire pyrolysis oil (TPO) obtained from the four runs of the pyrolysis plant were relatively similar to each other. TPO is a dark brown oil with a strong acrid smell. It interacts with human skin and it usually leaves spots and a bad smell for several days. So, careful treatment is required during use of TPO. The TPO properties for the four experiments are compared to commercial diesel and the results are tabulated in Table 3. It can be concluded that the density of TPO is more than the density of the commercial diesel and less than the density of the heavy fuel oil (985 kg/m<sup>3</sup> at 20 ◦C). The viscosity of the TPO is more than the viscosity of the commercial diesel and less than the viscosity of heavy fuel oil (200 cSt at 50 ◦C). Low viscosity of the tire pyrolysis oil is a wanted property in liquid to ease handling and transportation of the fuel. The flash temperature of a fuel gives an indication for the hazards of fire during storage and usage of the fuel. From Table 4, it can be concluded that the flash point of TPO is between 30 and 33 ◦C. The flash point of TPO is low if it is compared to conventional diesel fuel. The flash point for the conventional diesel fuel is 55 ◦C, for kerosene the flash point is 23 ◦C and for light fuel oil it is 79 ◦C. The flash point of the TPO is low because TPO is a mixture of un-refined liquids that has a wide range of distillation. Table 4 shows that the pour point of the TPO is relatively low when compared with diesel. The pH value of the TPO is between 4.27 and 4.37, which is considered low acidic fuel. This weak acidity is not a problem regarding storage and handling of the TPO, for example, soft drink companies like Cola or Pepsi utilize Polyethylene Terephthalate

(PET) plastic for storage and handling, although the value of pH for their drinks is about 2.5 [41]. Subsequently, TPO should be filtered, centrifuged and desulfurized before it can be utilized as a heavy furnace oil or blended with conventional diesel for furnaces, boilers and power plants.


**Table 4.** Characteristics of the tire pyrolysis oil (TPO) in comparison to commercial diesel.

ND: Not Determined, <sup>a</sup> at 40 ◦C, <sup>b</sup> at 30 ◦C.

#### *4.4. Characterization of the Pyrolysis Gas (Pyro-Gas) Fraction*

The compositional analysis showed that the pyrolysis gas (Pyro-Gas) consists mainly of light hydrocarbons and the concentrations of methane and hydrogen were the highest. Hydrogen concentration could reach 40% in volume at high pyrolysis temperatures of 550 ◦C [40]. Pyro-Gas has a high heating value (38–41 MJ/m3), which is considered close to the HHV of natural gas (around 39 MJ/m3), that could meet the need of energy during the pyrolysis process or other industrial applications. However, it would be very important to clean Pyro-Gas before burning to take out the hydrogen sulfide (H2S) from the Pyro-Gas fraction. The corrosive nature of H2S would damage the downstream equipment. H2S amount must be less than 500 ppm in order to be used for commercial alternative engines. In addition, SO2 emitted after the combustion of the gas fraction should satisfy the environmental constrains.

#### **5. Summary**

The gravimetric analysis of tire waste pyrolysis products from the reactor working at the optimum conditions is tire pyrolysis oil (TPO): 45%, pyrolysis carbon black (PCB): 35%, pyrolysis gas (Pyro-Gas): 10% and steel wire: 10%. It can be concluded that chlorine and sulfur were concentrated in the PCB fraction due to the nature of the tires' composition. Approximately all chlorine in the tires was concentrated in the PCB fraction. Therefore, it is very important to take this result into account once PCB is used as a fuel.

The main disadvantage of using the TPO as fuel is its strong acrid smell, which could last for a few days if it interacts with human skin or clothes. The density of TPO is more than that of conventional diesel. However, it is still less than the density of heavy fuel oil. On the other hand, the flash point of TPO is low if it is compared with the other refined fuels, which are considered a disadvantage associated with the hazards of fire during the usage and storage of the fuel.

The pH value of the TPO is considered to be of a weak acidic nature but this weak acidity is not a problem regarding storage and handling of TPO. The compositional analysis showed that the pyrolysis gas (Pyro-Gas) consists mainly of light hydrocarbons and the concentration of methane and hydrogen were the highest. The concentration of H2 is about 40% in volume and it has a High Calorific Value (38–41 MJ/m3), which could meet the process requirement of energy. However, it would be very important to clean Pyro-Gas before burning so as to take out the hydrogen sulfide (H2S) from Pyro-Gas, and hence reduce the acid rain problem.

The heavy fuel oil price in Palestine is around 0.5 \$/liter. However, the production cost of the tire pyrolysis oil (TPO) is around 0.15 \$/liter. This cost is less than one third the price of the conventional heavy fuel oil in Palestine. If three pyrolysis plants, such as the one in the present study with a capacity of 330 (tons/month), were founded, they would be able to get rid of the waste tires generated yearly in Palestine (3800 metric tons). If all the waste tires generated in Palestine were pyrolyzed, the annual import heavy fuel oil would decrease by 1700 tons (11,850 barrels) and produce more than 1330 tons of pyrolysis black carbon and 380 tons of steel wire, yearly. In addition, and most importantly, the environmental hazardous problem of the waste tires could be solved in the correct way, reducing the dependence on imported furnace oil and creating new job opportunities.

#### **6. Conclusions and Recommendations**

The following recommendations should be applied for more comfortable operation and safer work and environmental conditions:


**Author Contributions:** R.A., M.A. and A.J. conceived and performed the experiments used in the article; A.J., M.A. and R.A. analyzed the data; R.A., A.J., and M.A. wrote the paper. F.M.-A. and T.S. revised and supervised the manuscript. The authors shared the work on the structure and aims of the manuscript, paper drafting, editing and reviewing the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by An-Najah National University with grant number [ANNU-1819-Sc023].

**Acknowledgments:** The authors would like to acknowledge the Deanship of Scientific Research at An-Najah National University for their grant which helped in accomplishing the current study.

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

#### **List of Abbreviations**


#### **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* **Estimating the Optimum Tilt Angles for South-Facing Surfaces in Palestine**

#### **Ramez Abdallah 1, Adel Juaidi 1,\*, Salameh Abdel-Fattah <sup>1</sup> and Francisco Manzano-Agugliaro <sup>2</sup>**


Received: 7 December 2019; Accepted: 27 January 2020; Published: 1 February 2020

**Abstract:** The optimum tilt angle of solar panels or collectors is crucial when determining parameters that affect the performance of those panels. A mathematical model is used for determining the optimum tilt angle and for calculating the solar radiation on a south-facing surface on a daily, monthly, seasonal, semi-annual, and annual basis. Photovoltaic Geographical Information System (PVGIS) and Photovoltaic Software (PVWatts) is developed by the NREL (US National Renewable Energy Laboratory) are also used to calculate the optimum monthly, seasonal, semi-annual, and annual tilt angles and to compare these results with the results obtained from the mathematical model. The results are very similar. PVGIS and PVWatts are used to estimate the solar radiation on south-facing surfaces with different tilt angles. A case study of a mono-crystalline module with 5 kWP of peak power is used to find out the amount of increased energy (gains) obtained by adjusting the Photovoltaic (PV) tilt angles based on yearly, semi-annual, seasonal, and monthly tilt angles. The results show that monthly adjustments of the solar panels in the main Palestinian cities can generate about 17% more solar energy than the case of solar panels fixed on a horizontal surface. Seasonal and semi-annual adjustments can generate about 15% more energy (i.e., it is worth changing the solar panels 12 times a year (monthly) or at least 2 times a year (semi-annually). The yearly optimum tilt angle for most Palestinian cities is about 29◦, which yields an increase of about 10% energy gain compared to a solar panel fixed on a horizontal surface.

**Keywords:** optimal tilt angle; PV system; solar photovoltaic; solar radiation; PVGIS; PVWatts; Palestine

#### **1. Introduction**

Solar energy plays a crucial role in the future of sustainable energy as an environmentally friendly source of energy. Solar panel use has gained acceptance nowadays to become part of our modern-day life. The situation of energy in Palestine is different from that of surrounding countries due to several reasons: a scarcity of natural resources, the political conditions in Palestine are unstable, the difficult financial situation, and the high density of the population [1]. Palestine's natural resources and mineral wealth are scarce; this scarcity includes traditional sources of energy, such as oil and gas deposits. The high prices of hydrocarbon fuels in Palestinian cities are equivalent to those prices found in the most expensive cities in the world [1]. In addition, Palestine imports 100% of its fossil fuel and 87% of its electricity from other countries [1].

Palestine desperately needs all kinds of energy for growing its weak economy. There is no access to electricity for all Palestinian people for the whole day [1]. Palestinians strongly seek to overcome this dilemma by finding alternative sources of conventional fuel [2]. As in the case of the civilized world, they have found alternative sources of energy through which they hope to reduce their dependence on other countries. An example of this is the fuel crisis that the Gaza Strip experienced in mid-2013, which negatively affected the lives of people, especially in the health and transportation sectors; up to this time, Gaza Strip populations were living outside modern world era using primitive methods of mobility and lighting that disappeared a long time ago [2,3].

There is no doubt that the Palestinian tendency has shifted toward reliance on renewable sources of energy, which is consistent with the increasing global trends in the exploitation of alternative sources of energy. The trend in Palestine is concentrated in solar, biomass, wind, and geothermal energy [4]. The global trend is to find a clean, renewable, domestic sustainable energy source. This trend philosophy is increasing day by day, after significant damage caused by fossil fuels and their obvious risks to human health and the environment; toxic fumes and harmful gases from fossil fuel burning are a lingering health hazard endangering different forms of life on this planet. The most prominent environmental impact manifestations can be seen and felt right now, such as global warming, climate variability, ozone hole expansion, and acid rain, well as the frequent pollution of the seas by oil spills. The surge of oil prices is another driving force to find a replacement for conventional fossil fuels.

The Palestinians have succeeded to some extent in the exploitation of solar energy; they been using solar water heaters to heat water for domestic household uses and for limited industrial applications since the mid-seventies of the last century. The solar heater is a key component of every Palestinian home [5]. On the other hand, solar electricity generation remains limited due to research issues or limited activities of donors to help the population in disadvantaged areas. There are initial attempts at the exploitation of wind energy, as well as thermal energy. All these applications are still limited in the embryonic stage, with a promise of good potential [6,7].

#### *1.1. Physical Characteristics of Palestine*

Palestine, as it exists today, is formed from two separated land areas. The first is the West Bank with a total area of 5661 km2, and the other is the Gaza Strip with a total area of 362 km2 [8]. The West Bank is bordered to the north, west, and south by Israel and to the east by the Jordan River. It consists of eleven governorates: Jerusalem (East Jerusalem), Ramallah, Al-Bireh, Hebron, Nablus, Jenin, Tulkarm, Jericho, Tubas, Qalqiliya, Bethlehem, and Salfit. The Gaza Strip is located in the eastern Mediterranean Sea [8]. The Gaza Strip is bordered to the north and east by Israel, to the south by Egypt, and to the west by the Mediterranean Sea. It consists of five governorates: Rafah, Khan Yunis, North Gaza, Deir al Balah, and Gaza (as shown in the map in Figure 1).

**Climate:** Palestine has a Mediterranean climatic zone which has hot, long, dry summers and cool, short, wet winters [8]. However, the south of the Jordan Valley has a different climate zone and lies between the extreme desert conditions and dry steppe [8]. The climate within Palestine is affected by orographic effects, distance from the Mediterranean, altitude, and latitudinal position. Therefore, the climate of the West Bank can be divided into four main climatic zones: the Eastern Slopes, the Western Slopes, the Jordan Valley, and the Central Highlands. The climate of the West Bank, especially in the south, is affected by the vast nearby Arabian and the Negev deserts, especially during early summer and spring. Desert storms with hot winds full of dust and sand (locally known as Khamaseen winds) lead to increasing air temperatures and decreased humidity [8].

**Temperature:** The temperature of Palestine is comparatively high. The Jordan Valley and Jericho have the highest temperatures. The temperatures of the Jordan Valley are related inversely to altitude and decrease from south to north, with the highest temperature in the southern part at the Dead Sea. The Dead Sea is considered the lowest point on earth and is surrounded by a series of high mountains from both west and east, which create a natural greenhouse climate [8]. In summer, June, July, and August, the monthly average temperatures for the West Bank range from 20.8 ◦C to 30 ◦C. During winter months, December, January, and February, the monthly average temperatures in the West Bank range between 8.7 ◦C and 14.7 ◦C. The Gaza Strip lies between the semi-humid Mediterranean climate along the coast and the arid desert climate of the Sinai Peninsula. The monthly mean temperatures range from 13 ◦C in winter to 25 ◦C in summer [8].

**Sunshine Duration:** In Palestine, solar radiation varies from one city to another. June has the longest duration of sunlight, and December has the shortest duration of sunlight [8]. During the summer months, clear skies strengthen solar radiation, whereas during the winter, cloud cover reduces solar radiation [8].

#### *1.2. The Situation of Electric Energy in Palestine and Net-Metering Regulations*

The total electricity consumption for Palestine is 5137 GWh/year. About 87% is imported from the Israel Electric Corporation (IEC), while less than 8% is from the Palestine Electric Company. A small portion is imported from Jordan to supply Jericho in West Bank and from Egypt to supply Rafah in the Gaza Strip [1]. According to the Palestinian Central Bureau of Statistics (PCBS), more than 99.7% of the Palestinians in the West Bank and Gaza Strip are connected to the grid. For that reason, most of the solar PV in Palestine is an on-grid solar system [9].

The Palestinian Energy Authority (PEA) policy is to encourage the Palestinian people to utilize solar PV power to contribute to independency of Palestinian energy and support the Palestinian economy. To encourage solar PV energy, the PEA introduced the net-metering regulations.

The electricity generated by PV is supplied to the houses, where the grid serves as 'electricity storage'. In the case where solar PV produces more energy than that consumed by the houses, then the excess electrical energy is transferred to the grid. On the other hand, when the house electricity consumption is higher than the production, the house uses electricity from the grid. If the injection energy is more than the consumption during the month, 75% of the excess energy is credited to the next month. The credits are valid for 36 months. In times of higher consumption than production, the balance has to be paid by the consumer [10].

#### **2. Methodology**

Palestine is one of the countries enclosed between latitudes 40 degrees north and 40 degrees south (the so-called sun-belt). Palestine has more than 300 sunny days per year at an average brightness of 8 hours per day, increasing in summer and decreasing in winter. Previous studies show that the average daily solar radiation in Palestine reaches 5.4 kWh/m2, which is equivalent to an annual production of 1950 kWh of energy, which has enabled Palestine to be one of the best areas using solar energy, and the possibility of investing in it is economically feasible [5].

When installing photovoltaic cells in an area, one must first calculate solar radiation in this region to optimize installation layouts to maximize the generated electricity, and to calculate solar radiation first, we must recognize the so-called solar angles.

#### *Solar Angles*

**Declination angle (**δ**):** is the angle between the equatorial plane and the line that connects between the earth's center and the sun's center. The declination angle varies seasonally between 23.45 and −23.45 and is calculated from the following relationship.

$$\delta = -23.45 \cos\left( (n + 10.5) \frac{360}{365} \right) \tag{1}$$

n = the day number, such that n = 1 on the 1st of January.

**Latitude angle (**∅**):** is the angle between the equator and the line passing through the earth's center and a point on the earth's surface; northern latitudes are positive, and southern latitudes are negative.

**Hour angle (h):** is the angle between the plane containing the axis of the earth and the zenith and the other plane containing the axis of the earth and a point on the Earth's surface.

**The altitude angle (**α**):** is the angle between the horizontal and the rays of the sun. Therefore, the altitude angle for sunrise and sunset equals zero. The altitude angle (α) is calculated by [11]:

$$
\sin \alpha = \sin \delta \sin \mathcal{Q} + \cos \delta \cos \mathbf{h} \cos \mathcal{Q}.\tag{2}
$$

Note that the hour angles at sunrise and sunset (hs) have the same magnitude and can be calculated by substituting α = 0 into the previous equation (Equation (2)) to obtained:

$$\mathbf{h}\_{\mathbb{A}} = \cos^{-1}[-\tan(\mathcal{Q})\tan(\delta)].\tag{3}$$

**The tilt angle (**β**):** is the angle between the horizontal and the tilted surface.

The total solar radiation is the sum of the direct, diffuse, and reflected solar radiation. Firstly, direct radiation depends on the atmosphere's thickness, the quantity of water vapor, and the level of pollution in the atmosphere. Secondly, diffuse solar radiation is defined as all the solar radiation that comes from the sky except for the part that comes directly from the sun. Therefore, it is difficult to calculate because it totally depends on highly variable cloudiness and atmospheric clarity. Finally, reflected radiation from buildings, cars, the ground, and streets must be calculated individually. For that reason, reflected radiation is usually neglected. In this paper, an "ASHRAE standard atmosphere" is assumed. This assumption will overestimate the solar radiation. On the other hand, this paper

studies the optimum tilt angle rather than the available radiation. The total radiation obtained by the "ASHRAE standard atmosphere" assumption will be more than the actual available radiation, but the optimum tilt angles to achieve maximum radiation from any tilted surface will be the same. At the end of this study, the results will be compared to previously published papers.

In the Northern Hemisphere, it is considered that optimum orientation is due south [12,13]. Moreover, in 2013, Jafarkazemi and Saadabadi [14] found that the optimum orientation for Abu Dhabi was due south. Also, in 2010, Pavlovi´c, et al. [15] estimated the optimum orientation for Serbia to be due south. Finally, in 2018, Omar and Mahmoud [16] confirmed that the optimum orientation for Palestine is due south. Based on the previous studies, the present mathematical model considers the titled surfaces to be due south.

The daily total extraterrestrial radiation on a tilted surface with an angle surface β and facing south is [9]:

$$I\_d = \frac{24}{\pi} I\_o \left[ 1 + 0.034 \cos(\frac{2\pi\mathbf{n}}{365}) \right] \times \left[ \cos(\mathcal{Q} - \beta) \cos(\delta) \sin(\mathbf{h}\_s) + \mathbf{h}\_s \sin(\mathcal{Q} - \beta) \sin(\delta) \right]. \tag{4}$$

From Equation (4), for a particular day n to have an optimum tilt angle (β*op*,*d*), the total radiation Id must be derived with respect to the tilt angle (β) and equalized with the resulting derivative set to zero, i.e., *dId <sup>d</sup>*<sup>β</sup> = 0.0. to obtain:

$$\beta\_{op,d} = \left. \oslash - \tan^{-1} \right[ \frac{h \text{s}}{\sin(h \text{s})} \tan(\delta)]. \tag{5}$$

It could be concluded that adjusting the surface tilt angle every day is not practical, but it can be changed monthly or seasonally; the total radiation for a particular period can be obtained as:

$$I\_p = \sum\_{n=n\_1}^{n=n\_2} I\_{d\_{n\_1}} \tag{6}$$

where p is the particular period, n1 and n2 are the first and last days of a particular period, respectively. For obtaining the optimum tilt angle β*op*,*<sup>p</sup>* for a particular period (from n1 to n2), Ip must be derived with respect to the tilt angle (β) and equalized with the resulting derivative set to zero, i.e., *dIp <sup>d</sup>*<sup>β</sup> = 0.0 to obtain:

$$\beta\_{op,p} = \left. \oslash - \tan^{-1} \right[ \frac{\sum\_{n=n\_1}^{n=n\_2} \frac{24}{\pi} I\_o \left[ 1 + 0.034 \cos \left( \frac{2\pi n}{365} \right) \right] \sin(\delta) \text{hs}}{\sum\_{n=n\_1}^{n=n\_2} \frac{24}{\pi} I\_o \left[ 1 + 0.034 \cos \left( \frac{2\pi n}{365} \right) \right] \cos(\delta) \sin(\text{hs})} \text{.} \tag{7}$$

**A MATLAB program** (An-Najah National University Computer Labs.) was created to find the optimum tilt angle for any particular period using Equations (1)–(7). Using this program, the optimum tilt angle and total extraterrestrial radiation can be obtained for daily, monthly, seasonal, and yearly values. The major solar angles are shown in Figure 2a,b.

A comparison of the optimum tilt angle of the previous mathematical model was performed with the optimum tilt angle derived from the National Renewable Energy Laboratory (NREL) PVWatts program (NREL, 2017) [17] and PVGIS has been developed at the European Commission Joint Research Centre, at the JRC site in Ispra, Italy since 2001 [18]. The PVWatts Calculator is an online calculator that gives the amount of global radiation and the energy production of the PV using solar resource data for fixed locations throughout the world. The PVGIS gives the amount of global radiation and the energy production of the PV using the solar radiation data produced from satellite images. For the version of PVGIS used in this study, the satellite data used for the solar radiation estimates are from the METEOSAT satellites covering Europe, Africa, and most of Asia. Depending on satellite images, these images are captured every 15 minutes or 30 minutes [18]. Utilizing PVGIS, one image per hour

was used; this means that PVGIS can be used for any location regardless of the distance to a particular meteorological station.

**Figure 2.** Major solar angles: (**a**) Tilt angle (β) and Altitude angle (α) (**b**)Declination angle (δ), Latitude angle (∅), and Hour angle (h).

PVWatts and PVGIS are used to calculate the solar radiation on a horizontal surface with several tilt angles. A case study of 5 kWP of peak power is used to demonstrate the importance of the optimum tilt angles.

#### **3. Literature Review**

With a constant orientation facing south, one of the most important factors that affects the performance of a solar panel or collector is its tilt angle since changing the tilt angle changes the amount of solar radiation received by the solar panel or collector.

The performance of solar equipment (i.e., solar collectors for water and air heating, power generation, photovoltaic systems, etc.), is closely related to the inclined angle from the surface that absorbs the light (i.e., with solar equipment placement plane) [19].

The maximum daily energy can be achieved using solar tracking systems. A mechanical device should be added to follow the direction of the sun on its daily sweep. These devices are expensive, need additional energy for operation, and are not always applicable [20,21]; therefore, it is often more practical to direct the surface at an optimum tilt angle, β. Moreover, according to the Palestinian Central Bureau of Statistics, there are no solar tracking projects in Palestine [9].

The optimum tilt angle must be changed daily, monthly, seasonally, semi-annually or annually to maximize the amount of solar radiation received by solar equipment [22].

Solar energy should be given more attention because it its environmentally friendly with no pollution and is a sustainable source of energy. The amount of solar radiation changes with the latitude of the location, the day number through the years, and time of a day. There are still lots of errors during the installation of photovoltaic (PV) panels; significant errors are made by companies and personnel. These errors are wrong PV cell selection, wrong installation spaces and orientation of the PV panels, and also wrong tilt angles [23]. Jiménez-Torres, et al. stated that information of the available solar resource in a particular location on a surface in any position (inclination and orientation) is essential [24]. Mousavi, et al. mentioned that the efficiency of a PV panel can be optimized if the panel is positioned in such a way that it receives the maximum amount of incident solar radiation [25]. Table 1 provides the optimal tilt angle results in different countries of the world.


**Table 1.** Optimal tilt angle results in different regions of the world.


**Table 1.** *Cont*.

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

Since more than 99.7% of the Palestinian people in the West Bank and Gaza Strip are connected to the grid and the majority of the PVs in Palestine are an on-grid solar system in which the grid serves as electricity storage [9], the optimum tilt angle is important to produce the maximum energy generation during the month regardless of the peak load time during the day. The advantages of on-grid PV solar systems include the benefits of getting rid of the difficulties and energy losses associated with energy storage and the peak load consumption vs. the maximum power generation from the PV system.

Considering that most Palestinian communities are connected to the grid and the electricity tariff is constant during the day, the maximum energy production during the year is a vital factor for PV applications. In the present study, the optimum tilt angle is calculated to achieve this maximum energy production.

The daily, monthly, seasonal, semi-annual, and annual optimum tilt angles are calculated using the mathematical model. Also, a monthly, seasonal, semi-annual, and annual optimum tilt angle is calculated for two Palestinian cities, Jerusalem in the West Bank and Gaza city in the Gaza Strip, using PVWatts and PVGIS.

According to PEA regulations, 5 kWP is the maximum allowed peak power PV for houses connected to the grid in Palestine. The annual energy consumption for most residential houses does not exceed 5 kWP PV. Therefore, most houses implement 5 kWP. Moreover, there are no specific standards or specifications for selecting PV modules, but all the PV modules in Palestine are a polycrystalline or monocrystalline type. For that reason, in the present study, 5 kWP will be taken as a case study.

#### *4.1. Part One: Mathematical Model Results*

#### 4.1.1. Daily Optimum Tilt Angle

Equations (1), (3), (4), and (5) are used in the MATLAB program to calculate the daily optimum tilt angle βopt,d and the daily total extraterrestrial radiation for different latitudes and the main Palestinian cities. The results obtained from the MATLAB program are tabulated for recommended average days of months in Table 2 for different latitude angles and in Table 3 for the main Palestinian cities and are plotted in Figure 3 for different latitude angles as well as in Figure 4 for Jerusalem and Gaza city. The results show that the daily optimum tilt angle βopt,d increases with the increasing latitude of the selected location, and it changes with the number of day which takes the form of a triangular cosine function. Moreover, the optimum daily angle βopt,d for 21 March is slightly more than the latitude of the location whereas βopt,d for 22 March is slightly less than the latitude of the location. On the other hand, the optimum daily angle βopt,d for 20 September is slightly less than the latitude of the location whereas βopt,d for 21 September is slightly more than the latitude of the location. The same results were obtained for the main Palestinian cities. Daily optimum tilt angles βopt,d for the Palestinian cities are plotted in Figure 3. It can be concluded that negative values for βopt,d for the main Palestinian cities were obtained for the time period from n = 134 to 210.

**Figure 3.** The daily optimum tilt angle βopt,d for different latitude angles.

#### 4.1.2. Monthly Optimum Tilt Angle

Equations (1), (3), (4), (6), and (7) are used in the MATLAB program to calculate the monthly optimum tilt angle (βopt,m) and the monthly total extraterrestrial radiation (Im) for different latitude angles and for the main Palestinian cities. The results obtained from the MATLAB program are tabulated in Table 4 for different latitude angles and are plotted in Figure 5. It important to mention that when the optimum monthly tilt angle βopt,m is negative for a particular month, then the monthly total extraterrestrial radiation for that month is recalculated by considering that the optimum monthly tilt angle βopt,m is equal to zero. The reason for that is to determine if it is worth designing a PV or solar collector with a negative tilt angle. The negative sign means that the solar panel or collector is north facing, while the positive sign means that the solar panel or collector is south facing. From Figure 4, it can be noted that the variation of the monthly optimum tilt angle can be described by a triangular cosine function. The results obtained for the main Palestinian cities are tabulated in Table 5. The βopt,m for Jerusalem is plotted in Figure 6. From Table 5 it can be concluded that βopt,m increases with increasing latitude angle. Moreover, for latitude 0, βopt,m has a negative value from April to September, while for latitudes 5, 10, and 15, βopt,m has a negative value from April to August. Also, for latitude 20, βopt,m has a negative value for May to August. Additionally, for latitudes 25 and 30, βopt,m has a negative value for May to July. For latitude 35, βopt,m has a negative value for June and July, and finally, for latitude 40, βopt,m has a negative value for only the month of June. The optimum tilt angle in March is slightly higher than the latitude but is slightly lower than the latitude in September.



**Table 3.** The daily optimum tilt angle βopt,d for recommended average days of months for the main Palestinian cities.


**Figure 4.** The daily optimum tilt angle βopt,d for Jerusalem (**a**) and Gaza City (**b**).

From Table 5, it can be concluded that for the main cities of Palestine, the monthly optimum tilt angle is between -7.0 (in June) to 62.0 (in December). Moreover, it can be concluded that for the main cities of Palestine, the optimum tilt angle for March is slightly higher than the latitude of the city (βopt,m = 35), whereas for September it is slightly less than the latitude of the city (βopt,m = 29◦). The yearly extraterrestrial radiation is calculated with β = 0 and β = βopt on a daily, monthly, seasonal, semi-annual, and annual basis. The tilt factor is defined as the ratio between the radiation on a tilted surface and the radiation on a horizontal surface in the same period. The tilt factor is calculated on a daily, monthly, seasonal, semi-annual, and annual basis, and the results are shown in Table 6. Moreover, the tilt factors for daily, monthly, seasonal, semi-annual, and annual optimum tilt angles for the main Palestinian cities are presented in Table 7. It can be concluded from Tables 6 and 7 that the tilt factor on a daily basis is almost the same as that on a monthly basis. This means that it worth changing the tilt angle monthly and there is no need to change it daily.





*Energies* **2020** , *13*, 623




I: The yearly extraterrestrial radiation in kWh/m2; F: The tilt factor.



*Energies* **2020** , *13*, 623

 1.7613

 3592.91

 1.7552

 3501.03

 1.7103

 3478.27

 1.6992

 3148.68

 1.5382

**Figure 5.** The monthly optimum tilt angle βopt,m for different latitudes.

**Figure 6.** The monthly optimum tilt angle βopt,m for Jerusalem, Palestine.

#### 4.1.3. Seasonal Optimum Tilt Angle

It is more practical to change the tilt angle quarterly rather than daily or monthly. Therefore, the tilt angle is changed once every three months. The Palestinian climate can be divided into four seasons: winter (December to February), spring (March to May), summer (June to August), and autumn (September to November). The seasonal optimum tilt angle and the seasonal total extraterrestrial radiation for different latitudes and for the main Palestinian cities are calculated. The seasonal variation

of optimum angle βopt,s, and the seasonal total extraterrestrial radiation Is for different latitude angles and for the main Palestinian cities are listed in Tables 8 and 9, respectively. Moreover, the seasonal variation of optimum angle βopt,s for different latitude angles and for Jerusalem are plotted in Figures 7 and 8, respectively. From Table 7, it can be shown that for most of the Palestinian cities, the solar radiation on a surface tilted by βopt,s is 25% more than the radiation falling on a horizontal surface (β = 0). Moreover, there is a decrease of less than 2% in the solar radiation on a surface tilted by the seasonal optimum tilt angle compared to a surface tilted by monthly optimum tilt angles for all the main Palestinian cities.

From Table 9, it could be concluded that the minimum and maximum solar radiation at a seasonally optimum tilt angle occurs in autumn and summer, respectively. For Jerusalem, the minimum and maximum radiation values are (947.13 kW h/m2) and (1016.15 kW h/m2), respectively.

**Figure 7.** The seasonally optimum tilt angle βopt,s for different latitudes.

**Figure 8.** The seasonally optimum tilt angle βopt,s for Jerusalem, Palestine.



**Table 9.** The seasonally optimum tilt angles βopt,s and extraterrestrial radiation Is in kWh/m2 for the main Palestinian cities.


#### 4.1.4. Semi-annually Optimum Tilt Angle

As mention above, some systems require a high efficiency of the solar system during winter. On the other hand, some systems require a high efficiency of the solar system during summer. Therefore, it is worth finding the semi-annual optimum tilt angle during the warm months and the cold months. The tilt angle is changed once in a period of six months. The warm period is from April to September while the cold period is from October to March. The semi-annual variation of optimum angle βopt,sa and the semi-annual total extraterrestrial radiation Isa for different latitude angles and for the main Palestinian cities are listed in Tables 10 and 11.

The results presented in Tables 6 and 7 show that the total yearly radiation received by a surface tilted by βopt,s on a seasonal basis is almost the same as that received by a surface βopt,sa on a semi-annual basis.

#### 4.1.5. Annually Optimum Tilt Angle

The yearly optimum tilt angles were calculated, and the yearly extraterrestrial radiation Iy gained by a surface tilted by the yearly optimum angle was calculated for different latitudes. The results are listed in Table 12. Similarly, for the main Palestinian cities, the results are listed in Table 13. Table 6 shows the tilt factor on a yearly basis. It could be noted that the increased solar energy on a surface tilted by a yearly optimum angle compared to the horizontal surface increases with increasing latitude of the location. The overall yearly increase of solar radiation is about 1.0% for the latitude of 10.0◦, increasing up to 53% for the latitude of 60.0◦.

From Table 7, it can be concluded that the increase in solar radiation for a surface tilted by a daily optimum angle relative to a surface tilted by a fixed yearly optimum angle is about 10.6% for most of the Palestinian cities, while the increase in solar radiation for a surface tilted by a monthly optimum angle is about 10.4%. This means that changing the angle each month gives approximately the same increase in solar radiation by changing it daily. Moreover, the increase in solar radiation for a surface tilted by a seasonally optimum angle relative to a surface tilted by a fixed yearly optimum angle is about 8.6% for most of the Palestinian cities, while the increase in solar radiation for a surface tilted by a semi-annually optimum tilt angle is about 8.4%. This means that changing the angle semi-annually gives approximately the same increase in solar radiation by changing it seasonally.

Finally, it can be shown that it is so important to tilt the collector or the panel with the optimum tilt angle. It is clear that for most Palestinian cities, the daily and monthly tilt angle energy gains reached 27% with respect to the horizontal surface, whereas the seasonal and semi-annual tilt angle energy gains reached 25%. Lastly, the energy gains for the yearly tilt angle were about 15% for most of the Palestinian cities.

#### *4.2. Part Two: Results Based on PVWatts and PVGIS*

Part one results showed a small variation between Palestinian cities by using either the optimum tilt angles as one approach or the amount of total radiation as the second approach. For that reason, two major cities in Palestine are selected; the first city is Jerusalem which is located in the West Bank and the other is Gaza city located in the Gaza Strip.

Jerusalem is located in the West Bank where there is a meteorological station; PVWatts depends on this station for global radiation. In Gaza city which is located in the Gaza Strip, the nearest meteorological station is Be'Er Sheva.

The reason for using PVWatts and PVGIS for this comparison is that PVWatts uses solar resource data collected from a ground station whereas PVGIS uses data produced from satellite images.

PVWatts does not estimate the optimum tilt angle, whereas PVGIS estimates only the yearly optimum tilt angle. Therefore, the optimum tilt angle for each city and particular time period is found by estimating PV outputs with various tilt angles for this particular period so the optimum tilt angle giving the maximum output is found.



**Table 11.** The semi-annual optimum tilts angle βopt,sa and extraterrestrial radiation Isa in kWh/m2 for the main Palestinian cities.





**Table 13.** The annual optimum tilt angle βopt,y and extraterrestrial radiation Iy in kWh/m2 for the main Palestinian cities.

The optimum tilt angles for two Palestinian cities, Jerusalem and Gaza city, are shown in Table 14 for the different time periods. Figures 9 and 10 show a comparison of monthly and seasonal optimum tilt angles for Jerusalem and Gaza using PVWatts and PVGIS with those calculated by the mathematical model. It is clear that the results are very close. This means that the mathematical model in part one gives accurate results for the optimum tilt angle.

**Table 14.** A comparison of monthly, seasonal, semi-annual, and yearly optimum tilt angles (degrees) for Jerusalem and Gaza City based on PVWatts and PVGIS with those calculated in the present work.


It worth mentioning that PVGIS and PVWatts do not estimate a negative tilt angle. This means the minimum tilt angle is zero (horizontal surface).

**Figure 9.** A comparison of monthly optimum tilt angles based on PVWatts and PVGIS with those calculated by the mathematical model: Jerusalem (**a**), Gaza City (**b**).

**Figure 10.** A comparison of seasonal optimum tilt angles based on PVWatts and PVGIS with those calculated by the mathematical model: Jerusalem (**a**), Gaza city (**b**).

#### **5. Case Studies of two Palestinian Cities**

To perform this study, the cities of Jerusalem and Gaza city have been selected. Jerusalem is located in the West Bank whereas Gaza city is located in the Gaza Strip.

PVGIS and PVWatts are used to calculate the energy generated in the two locations mentioned above for a photovoltaic composed of a monocrystalline silicon module with 5 kWP of peak power installed. The comparison is made with five configurations: a system with a PV in a horizontal position and annual, semi-annual, seasonal, and monthly optimum tilt angle positions. In this way, the differences and gains due to the optimal position can be quantitatively verified.

Firstly, PVGIS and PVWatts are used to obtain the monthly average global radiation for Jerusalem and Gaza. Tables 15 and 16 show the monthly average global radiation in Jerusalem and Gaza city for different tilt angles (β).


**Table 15.** The monthly average global radiation in Jerusalem for different tilt angles (β).

The maximum radiation value is reached in the horizontal PV in June for the two cities and using PVGIS and PVWatts, whereas the maximum monthly radiation for a PV tilted at an annual optimum tilt angle is achieved in July. It can be concluded that the annual average radiation (kWh/m2/day) for Gaza city is slightly higher than that for Jerusalem. Moreover, the annual average radiation (kWh/m2/day) calculated by using PVGIS is (6.12) slightly higher than the annual average radiation (5.96) calculated using PVWatts.


**Table 16.** The monthly average global radiation in Gaza city for different tilt angles (β).

It can be noticed that annual average radiation (kWh/m2/day) for a PV tilted at the yearly optimum tilt angle calculated by the mathematical model suggested in this study and the annual average radiation (kWh/m2/day) for a PV tilted at the yearly optimum tilt angle calculated by PVGIS and PVWatts is very similar. For example, the yearly optimum tilt angle calculated by the suggested mathematical model in this study for Jerusalem is 29◦ and that when using PVWatts is 26◦ as shown in Table 15. However, the annual average radiation is the same (5.85 kWh/m2/day) as shown in Table 16.

In another example, the yearly optimum tilt angle calculated by the suggested mathematical model in this study for Gaza is 29◦and that when using PVGIS is 28◦as shown in Table 16. However, the annual average radiation is the same (6.12 kWh/m2/day) as shown in Table 16.

For the city of Jerusalem and using PVGIS, it is observed that the annual average radiation on a horizontal plane is (5.55 kWh/m2/day) as shown in Table 15 while for 28◦(PVGIS yearly optimum tilt angle), it is (6.08 kWh/m2/day) as shown in Table 15, and for 29◦(yearly optimum tilt angle obtained from the mathematical model), it is (6.07 kWh/m2/day) as shown in Table 15. It can be concluded that increasing the tilt angle more than the optimal value also leads to a decrease in the annual average radiation, obtaining only (3.53 kWh/m2/day) as shown in Table 16 for a vertical plane β (90◦).

For the city of Gaza and using PVWatts, it is observed that the annual average radiation on a horizontal plane is (5.36 kWh/m2/day) as shown in Table 16 while for 29◦(the same yearly optimum tilt angle using PVWatts and the mathematical model), it is (5.96 kWh/m2/day) as shown in Table 16. For a vertical plane, it is only (3.5 kWh/m2/day) as shown in Table 16.

For a photovoltaic system with a peak power of (5 kWh) in Jerusalem, the following monthly values are calculated using PVGIS and PVWatts.

Tables 17 and 18 show the monthly energy generated in Jerusalem and Gaza city by a 5 kWh system.

Figures 11 and 12 show the monthly energy generated in Jerusalem and Gaza city by a 5 kWh system.

It can be concluded that the energy generated by a PV system fixed at a yearly optimum tilt angle is more than 10% higher than that for a horizontal PV system. Moreover, the energy generated by a PV system fixed at semi-annual tilt angle gives 15% more energy gain as shown in Table 19.


**Table 17.** Monthly energy generated by a 5 kWh system in Jerusalem.



**Figure 11.** Monthly energy generated in Jerusalem by a 5 kWh system.

**Figure 12.** Monthly energy generated in Gaza city by a 5 kWh system.

It can be noted that a slight energy gain can be obtained by a seasonally optimum tilt angle compared to a semi-annual optimum tilt angle.

Changing the tilt angle 12 times through the year to obtain the monthly optimum tilt leads to about 17% energy gain with respect to a horizontal PV system as shown in Tables 17 and 18. This means that it is worth fixing the PV system at the yearly optimum tilt angle and changing the optimum tilt angle at least two times per year.

#### **6. Comparison of the Present Work with Previously Published Papers**

Firstly, a comparison of the present mathematical model was made with H. Darhmaoui and D. Lahjouji [37]. They used a mathematical model to estimate the solar radiation on a tilted surface at 35 sites in different countries of the Mediterranean region. They found the yearly optimum tilt angle for Gaza city to be 32◦, whereas in the present mathematical model the value was 29◦.

Secondly, a comparison of the present mathematical model with the models of Nijegorodor et al. [39] and P. Talebizadeh et al. [19] was made. Nijegorodor et al. and P. Talebizadeh et al. suggested monthly and yearly equations for calculating monthly and yearly optimum tilt angles. Monthly and yearly equations were used to estimate the optimum monthly and yearly tilt angles for Jerusalem, and a comparison of these angles with those obtained by the present mathematical model is given in Table 19 and is plotted in Figure 13. Results show that the yearly optimum tilt angles are very close while the monthly optimum tilt angle is slightly different.


**Table 19.** A comparison of monthly and yearly tilt angles for Jerusalem based on Nijegorodor et al. and P. Talebizadeh et al. with the values calculated in the present mathematical model.

**Figure 13.** Comparison of monthly and yearly tilt angles for Jerusalem based on Nijegorodor et al [39]. and P. Talebizadeh et al [31], with values calculated in the present mathematical model.

#### **7. Conclusions**

It is clear that the suggested mathematical model gives accurate results for the optimum tilt angle compared to PVWatts and PVGIS. It is worth designing a solar panel or collector with a negative tilt angle for locations at latitudes from 0◦ to 30◦. For all the Palestinian cities, there is no need to design a collector with a negative tilt angle. Compared to a horizontal surface, yearly, semi-annual, seasonal, and monthly adjustments of solar panels in the main Palestinian cities will generate approximately 10%, 15%, 15%, and 17% more energy, respectively, whereas compared to the fixed optimum yearly tilt angle, monthly adjustments of the solar panels in the main Palestinian cities can generate about 6% more energy. Seasonal and semi-annual adjustments can generate about 5% more solar energy. This shows that solar panels or collectors have to be adjusted monthly or semi-annually. For the

roof top-solar PV rated power (5 kWP) which was taken as a case study for the present work, it is convenient to change it semi-annually. Moreover, this proves the importance of an accurate tilt angle for the solar panels or collectors. If a supporting structure for the PV is well designed, then the solar PV can be easily adjusted in a limited time and can be performed by an unqualified person.

The daily optimum tilt angle for the main Palestinian cities ranges from −8.0◦ in June to 62◦ in December. The monthly optimum tilt angle for the main Palestinian cities throughout the year is between −7.0◦ in June and 62◦ in December. The seasonal optimum tilt angle for the main Palestinian cities throughout the year is between a minimum value of −1.0◦ in June and a maximum value of 58◦ in December. Moreover, the optimum tilt angle for autumn is more than the latitude of the city (44◦) and for spring, less than the latitude of the city (17◦). The semi-annual optimum tilt angle for the main Palestinian cities throughout the year is between 7.0◦ in the warm period (April to September) and 52◦ in the cold period (October to March).

The yearly average radiation in Jerusalem using PVGIS is (5.55 kWh/m2/day) based on a horizontal surface while using a panel tilted with a yearly optimum tilt angle, the value is (6.08 kWh/m2/day). The yearly average radiation in Jerusalem using PVWatts is (5.37 kWh/m2/day) based on a horizontal surface while using a panel tilted with a yearly optimum tilt angle, the value is (5.85 kWh/m2/day). The yearly average radiation in Gaza using PVGIS is (5.57 kWh/m2/day) based on a horizontal surface while using a panel tilting with a yearly optimum tilt angle, the value is (6.12 kWh/m2/day). The yearly average radiation in Jerusalem using PVWatts is (5.36 kWh/m2/day) based on a horizontal surface while using a panel tilted with a yearly optimum tilt angle, the value is (5.96 kWh/m2/day). It is clear that the PVGIS results are slightly higher than the PVWatts results for Palestinian cities.

The yearly optimum tilt angle for most Palestinian cities is about 29◦. It can be concluded that the monthly optimum tilt angle for a particular month is approximately the average of the daily optimum tilt angle for the particular month, while the seasonal, semi-annual, and yearly optimum tilt angles are approximately the average monthly optimum tilt angle for that particular period.

**Author Contributions:** R.A., and A.J. conceived of the idea and wrote the article; A.J. and R.A. analyzed the data; R.A., A.J., and S.A.-F. wrote the paper. F.M.-A. supervised the research and revised the manuscript. All authors contributed to the structure and aims of the manuscript, paper drafting, editing, and review. All authors have read and approved the final manuscript.

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

**Acknowledgments:** The authors would like to acknowledge the Palestinian Central Bureau of Statistics (Eng. Haitham Zeidan) and An Najah National University for facilitating this research.

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

#### **Abbreviations**



#### **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* **A Rebalancing Strategy for the Imbalance Problem in Bike-Sharing Systems**

#### **Peiyu Yi, Feihu Huang and Jian Peng \***

College of Computer Science, Sichuan University, Chengdu 610065, Sichuan, China **\*** Correspondence: jianpeng@scu.edu.cn; Tel.: +86-28-8546-9720

Received: 25 May 2019; Accepted: 1 July 2019; Published: 4 July 2019

**Abstract:** Shared bikes have become popular traveling tools in our daily life. The successful operation of bike sharing systems (BSS) can greatly promote energy saving in a city. In BSS, stations becoming empty or full is the main cause of customers failing to rent or return bikes. Some truck-based rebalancing strategies are proposed to solve this problem. However, there are still challenges around the relocation of bikes. The truck operating costs also need to be considered. In this paper, we propose a customer-oriented rebalancing strategy to solve this problem. In our strategy, two algorithms are proposed to ensure the whole system is balanced for as long as possible. The first algorithm calculates the optimal state of each station through the one-dimensional Random Walk Process with two absorption walls. Based on the derived optimal state of each station, the second algorithm recommends the station that has the largest difference between its current state and its optimal state to the customer. In addition, a simulation system of shared bikes based on the historical records of Bay Area Bikeshare is built to evaluate the performance of our proposed rebalancing strategy. The simulation results indicate that the proposed strategy is able to effectively decrease the imbalance in the system and increase the system's performance compared with the truck-based methods.

**Keywords:** bike sharing; energy saving; system rebalancing

#### **1. Introduction**

Energy saving is a vital task for the sustainable development of a city [1]. The government attaches great importance to the development of bike sharing systems (BSS), since they can greatly help reduce urban energy consumption [2,3]. We think the essential condition for the successful operation of BSS is solving the imbalance problem.

In BSS, the imbalance problem means some stations are empty (there are no available bikes for renting) while some stations are full (there are no empty docks for returning). A customer will experience an unsatisfactory service when there are no available bikes or docks at a station. Improving customer satisfaction is one of the tasks for BSS operators.

In order to deal with the imbalance problem, the operators of BSS need to relocate bikes among stations by trucks. Some static truck-based strategies for relocating bikes have been proposed [4–6]. The repositioning time is usually at the system's idle time (such as midnight) in these truck-based strategies. Recently, some work has paid attention to dynamic truck-based strategies which relocate bikes based on the advanced prediction of bikes [7–9]. However, the truck-based approaches are not only expensive to operate but also violate the low-carbon goal of bike-share due to the amount of pollution from the trucks.

There are also some strategies that encourage customers to participate in bike allocation. For example, some bike-sharing systems incentivize customers to rent bikes from stations with many available bikes and return to those with many empty docks by giving fare discounts [10–12]. Implementing the customer-oriented rebalancing strategies should consider the following tasks:

(1) Determine the optimal number of bikes for each station; (2) Deal with the changing of bike demand over time in a day; (3) Suggest the incentive station which cannot be too far from customer's original target station; and (4) Set an appropriate bonus for the customer under a limited budget. In existing studies, the incentive strategies will be implemented only when a station's number of bikes is not enough. However, the station is already in a bad state when its number of bikes is not enough, which will effect the performance of incentive measures.

In this paper, we propose a strategy which can help the system choose a suitable time and determine a suitable station for implementing the incentive measure. Based on the proposed strategy, a bike-sharing system can effectively lengthen its uptime. First, the proposed strategy calculates the most suitable number of available bikes for each station. Then it determines the station which has the largest capacity for renting. In addition, the walking distance is also considered in the strategy. Take an instance in Figure 1—a customer intends to rent a bike at station S1. At present, S1 has eight available bikes, which is less than the optimal state (ten available bikes). According to the idea of our proposed strategy, a returning action is helpful for S1 to reach its optimal state rather than a renting action. Therefore, the system should incentivize him to rent a bike at station S2. S2 is within r1 meters away from S1. After this rebalancing, the states of S1 and S3 will not be worse. At the same time, S2 can also reach its optimal state. For the bike-sharing system, the result reduces the probability of the imbalance problem occurring in future. The main contributions of this work are listed as follows.


**Figure 1.** Example of renting process.

The rest of this paper is organized as follows. Section 2 gives an outline of the development of BSS and research on rebalancing optimization; the problem statement and the framework are presented in Section 3; Section 4 introduces the details of the RWTAW and TLTF algorithms; we talk about the performance of the proposed rebalancing strategy in Section 5; finally, we conclude this work and look forward to future research in Section 6.

#### **2. Related Work**

BSS has achieved great development, with the goal of solving the "last kilometer" problem in cities [13], since it was first built in Amsterdam in 1965. The first kind of BSS is the traditional dock-based bike sharing system in which a customer must rent and return the bike at a dock and the number of docks at each station is fixed, such as the CitiBike in New York and the Hubway in Boston. The second is the dockless bike-sharing system such as the Mobike and Ofo in China. Each bike has a unique QR code and GPS tracking model. These two technologies increase the ease of use and management. The biggest advantage of the dockless bike-sharing system is customers can easily rent bikes via a smartphone app and park the bikes at any valid places [14–16]. The convenience also brings some problems, for example, operators have to rebalance the bikes. The existing studies on BSS mainly focus on the following aspects.

#### *2.1. System Rebalancing*

In BSS, system rebalancing is the main method for solving the imbalance problem. The truck-based rebalancing strategies and the customer-oriented rebalancing strategies are two popular approaches. In this section, we will introduce these two approaches in detail.

#### 2.1.1. Truck-Based Rebalancing

The truck-based rebalancing approaches have two main challenges, that is, determining the optimal inventory for each station and designing the optimal truck route with budget constraint. The approaches can be summarized into two categories according to the time of rebalancing. (1) Static Repositioning. Chemla et al. [4] used a branch-and-cut algorithm based on the tabu search to obtain the upper bounds and feasible routing solutions. However, their work only focused on solving the repositioning with a single vehicle in part of the city area. For the static bicycle repositioning problem under multiple vehicles, Raviv et al. [17] used the mixed integer linear program (MILP) to minimize the routing cost and customer dissatisfaction. Schuijbroek et al. [18] presented a rigorous cluster-first route-second heuristic to solve the inventory determining and routing decision problems. For the dockless system, Liu et al. [19] developed an enhanced version of chemical reaction optimization to solve the static bike repositioning problem by considering multiple heterogeneous vehicles, multiple depots and multiple visits. Pal and Zhang [5] presented a Novel Mixed Integer Linear Program for solving the static rebalancing problem in a free-floating bike sharing system. They also presented a hybrid nested large neighborhood search for large-scale bike sharing programs. (2) Dynamic Repositioning. Chiariotti et al. [20] used the Birth-Death Process to model a station's occupancy and proposed a dynamic rebalancing strategy. They first determined the optimal inventory for stations and then modeled the path-planning problem as a vehicle routing problem and adopted the greedy heuristic to dynamically redistribute bikes every hour. Contardo et al. [8] presented a mathematical programming formulation for the dynamic repositioning problem and obtained the lower bounds and feasible solutions by using the decomposition schemes. Li et al. [9] proposed a spatio-temporal reinforcement learning (STRL) model for dynamic bike repositioning. They first proposed an inter-independent inner-balance clustering algorithm to cluster stations into groups. Then they utilized the O-Model and the I-Model to predict the demands of renting and returning respectively. Lastly, for each cluster, the STRL obtains the optimal inner-cluster reposition policy based on a deep neural network. Liu et al. [21] developed a multi-source optimization approach for the rebalancing problem. They proposed an Adaptive Capacity Constrained K-centers Clustering

(AdaCCKC) algorithm which can reduce the large-scale multiple vehicle routing problem to an inner cluster one vehicle routing problem with guaranteed feasible solutions. Then the vehicle rebalancing route can be optimized by a mixed integer nonlinear programming (MINLP) model.

The drawbacks of the static and dynamic truck-based strategies are the expensive truck operating costs and the pollution of tail gas. Some researchers therefore now pay more attention to the customer-oriented rebalancing strategies.

#### 2.1.2. Customer-Oriented Rebalancing

The customer-oriented rebalancing approach is to solve the imbalance problem by providing customer incentives (bonus, fare discount or otherwise) without deploying trucks. For example, Singla et al. [11] proposed a method on the basis of a crowdsourcing mechanism [22] to incentivize customers to rent or return bikes at specific station. In their work, the dynamic pricing mechanism with a limited budget was proposed. It utilizes the regret minimization approach to maximize rebalancing efficiency. Chemla et al. [10] proved that the dynamic regulation problem is NP-hard (non-deterministic polynomial-time) and presented a heuristic to circumvent it. In their work, the pricing technique based on the linear programming was proposed. Fernandez et al. [23] designed a bike sharing system simulator to evaluate incentive-based rebalancing strategies. Some incentive strategies have been adopted in BSS. For instance, the work of Fricker and Gast [12] was used in the Velib+ system (in Paris). The incentive strategy is that customers will receive additional riding time for free if they return bikes at uphill stations. They also proved that the system can improve the operating efficiency if the customers who enjoy the incentive strategy park their bikes at two or more stations in a bad state. According to the current and predicted state of the system, Pfrommer et al. [24] proposed the price incentive scheme which takes into account the model-based predictive control principles. In a free-floating bike sharing system, Pan et al. [25] proposed a novel deep reinforcement learning framework for incentivizing users. Their work modeled the rebalancing problem as a Markov decision process and took both spatial and temporal features into consideration. The proposed Hierarchical Reinforcement Pricing (HRP) algorithm shows great performance on a dataset from Mobike, a major Chinese dockless bike sharing system. More recently, Diez et al. [26] proposed a persuasion model to recommend the optimal station and the shortest riding route between two locations. The model combines the argumentation theory and the characteristics of customers. In our work, the rebalancing strategy, which is based on the Random Walk process and integrates the walking distance between stations, was designed. The proposed strategy helps to lengthen the station's working time by inferring the optimal station for customers.

#### *2.2. System Prediction*

The prediction of bike demand is crucial for the bike repositioning problem. An accurate prediction can improve the effect of repositioning. Some station-level predicting models have been proposed. For example, Zeng et al. [27] proposed a station-centric model which takes into account the global features such as weather, user activity and season. Kaltenbrunner et al. [28] adopted the auto-regressive moving-average (ARMA) model to predict the number of bikes and docks for each station. Similarly, Yoon et al. [29] used a modified autoregressive integrated moving average (ARIMA) model to predict the available resources at each station. In their work, the spatial interaction and temporal factors are considered. Unlike the station-level prediction, some models are proposed according to the clustering of stations with similar patterns of usage. For example, Li et al. [30] proposed a hierarchical prediction model to predict the check-out/in of each station cluster. They first clustered the stations by the bipartite clustering algorithm based on the geographical locations of stations and transition patterns and then predicted the entire traffic in the city by way of the gradient boosting regression tree (GBRT) model. Finally, the rental proportion across clusters was predicted by a multi-similarity-based inference model. Chen et al. [31] proposed a dynamic cluster-based framework for over-demand bike prediction. Recently, neural networks have been adopted to predict bike demand. Liu et al. [32] proposed

a prediction model based on artificial neural networks (ANN). The model utilizes a set of features such as human mobility data, POI and station network structures. Lin et al. [33] proposed a novel graph convolutional neural network to predict the hourly demand in a BSS.

#### *2.3. System Design and Traffic Pattern Analyzing*

System design, including the station layout, capacity, and price policy, is the premise of successfully operating a bike sharing system. Dell'Olio et al. [34] calculated the potential demand of bike usage in a city and the price the customer would like to pay. They also proposed a location model with the help of geographical information about the stations. Chen and Sun [35] proposed a mathematical model to formulate the layout of public bike stations with the objective of minimizing users' total travel time. Lin and Yang [36] designed the strategic planning of BSS with service level considerations. They proposed an optimization model including a nonlinear integral and greedy heuristic to help operators determine the number of stations, the network structure of bike paths connecting the stations and the travel paths for users between each origin and destination station. In addition to the system design and efficient operation, the traffic pattern analysis is another important topic, such as the bike usage pattern and the transition pattern across stations in BSS [37–39]. Understanding the traffic pattern of BSS is helpful for operating the system efficiently and knowing the mobility of a city. Furthermore, the massive real trajectories generated by users help to solve some city problems, such as bike lanes planning [40] and illegal parking detection [41].

#### **3. Design Overview**

#### *3.1. Problem Statement*

In this paper, we focus on a bike sharing system in which the station is fixed. Generally, a BSS has *n* stations. The set of stations is defined as *S* = {*S*1, *S*2, ... , *Sn*}. The maximum capacity of station *Si* is defined as *Mi*. *mi*(*t*) represents the available bikes at station *Si* at time *t*. We model the bike sharing system as a fully connected graph *G* = {*S*, *I*}. Each node *Ii*,*<sup>j</sup>* = (*Si*, *Sj*) corresponds to the edge between two stations *Si* and *Sj*. The edges are weighed by walking distance metric *ri*,*j*. Notations used in this paper are shown in Table 1.

We set the scene of rebalancing as follows. A customer will query the information about stations on his mobile phone before his trip. After the customer selects the station at which he will rent a bike, the bike system immediately checks the station's state. If the number of the station's available bikes is less than its optimal state, the system would incentivize the customer to another station according to the states of nearby stations and the walking distance. We hope the whole system can obtain the least values of three metrics, that is, the unworking time, the proportion of lost customers and the times of reposition. The detailed definitions of the three metrics are introduced in Section 5. The selection of the recommended station at which the system incentivizes the customer to rent bikes is the main task in this work. Formally, the problem is defined as follows:

Problem Definition: Given the station *Si* at which a customer will rent a bike, the bike system will determine whether to incentivize the customer to another station *Sj* or not. If needed, the bike system should output the recommendation station *Sj*.


**Table 1.** Notations.

#### *3.2. Framework*

Figure 2 shows the framework of our bike-sharing system with the rebalancing strategy, which consists of the following three components:


**Figure 2.** Framework of simulate system.

#### **4. Methodology**

#### *4.1. The Optimal State Calculating*

In the BSS, we hope the station can keep the optimal state as long as possible, thus it should have the lowest possibility that it will be empty or full in future. Motivated by this idea, *OPi*,*<sup>k</sup>* is defined as the optimal state of *Si* at time slice *k*. It can be calculated on the basis of the probability that station becomes either empty state (*mi*(*t*) = 0) or full state (*mi*(*t*) = *Mi*). Formally, it can be defined as follows:

$$OP\_{i,k} = \underset{a}{\text{arg min}} \, P\_{a\prime} \, \forall a \in \{1, 2, \dots, M\_i - 1\} \tag{1}$$

where *Pa* is the probability that station becomes either empty or full during time *t*(*k*) and it is calculated by:

$$P\_a = P\_{a \to 0} + P\_{a \to M\_i} \tag{2}$$

where *Pa*→<sup>0</sup> and *Pa*→*Mi* respectively represent the probabilities that station *Si* becomes empty state and full state when the current state *mi*(*t*) = *a*.

Random Walk is a mathematical statistical model [44,45], which has a wide range of applications. For example, researchers have studied the application of the Random Walk model in finance prediction [46], high level data classification [47] and social network optimization [48,49]. A one-dimensional random walk can be looked at as a Markov chain whose state space is given by the integers *τ* (*τ* = 0, ±1, ±2, . . . ). The Random walk with two absorption states (*τ*<sup>1</sup> and *τ*2) at both ends is an important model. The model stops walking once it reaches either absorption state. For an initial state *τ<sup>a</sup>* (*τ*<sup>1</sup> < *τ<sup>a</sup>* < *τ*2), the probability *Pτ<sup>a</sup>* represents walking from *τ<sup>a</sup>* to either absorption state. *Pτ<sup>a</sup>* can be calculated as follows:

$$P\_{\tau\_d} = \frac{(\frac{q}{p})^{\tau\_d - \tau\_1} - (\frac{q}{p})^{\tau\_2}}{1 - (\frac{q}{p})^{\tau\_2}}\tag{3}$$

where *p* represents the probability of walking toward *τ*<sup>2</sup> at each step, *q* represents the probability of walking toward *τ*<sup>1</sup> at each step.

In this paper, we describe the state changing process of each station by the one-dimensional Random Walk process. In a bike sharing system, a station's state has only two changing directions: the empty state (*mi*(*t*) = 0) and the full state (*mi*(*t*) = *Mi*). The state of a station goes one step toward the empty state when a renting event occurs. In the same way, the state of station goes one step toward the full state when a returning event occurs. For a middle state *mi*(*t*) = *a*, if 0 < *a* < *OPi*,*k*, the smaller the value of *a*, the larger the probability that the station will become empty; if *OPi*,*<sup>k</sup>* < *a* < *Mi*, the larger the value of a, the larger the probability that the station will become full. We assume that the station will stop working and wait for supplementary as long as its state becomes either empty or full. Thus, the empty state and full state of a station correspond to two absorption walls.

Specifically, given the total changing steps *e*<sup>1</sup> and *e*2, if the station becomes empty after *e*<sup>1</sup> steps, the probability for this case is calculated by the following formula:

$$P\_{a \to 0} = \sum\_{i=0}^{\varepsilon\_1} \mathbb{C}\_{a+2i}^i \times P\_{rent}^{a+i} \times P\_{return}^i \tag{4}$$

Similarly, *Pa*→*Mi* is given as:

$$P\_{a \to M\_i} = \sum\_{i=0}^{c\_2} C\_{M\_i - a + 2i}^i \times P\_{return}^{M\_i - a + i} \times P\_{rent}^i \tag{5}$$

*Prent* and *Preturn*, respectively, represent the probabilities of a renting event and a returning event for each step at a station. The constraints of total steps *e*<sup>1</sup> and *e*<sup>2</sup> are:

$$e\_1 = \min\{\operatorname{Rect}\_{i,k} - a\_\prime \operatorname{Return}\_{i,k}\} \tag{6}$$

$$\epsilon\_2 = \min\{a - M\_l + \operatorname{Return}\_{i,k} \operatorname{Rent}\_{i,k}\}\tag{7}$$

In a bike sharing system, the process by which customers arrive at a station to rent or return bikes can be described by the Poisson process [43]. The renting intensity of Poisson process *Renti*,*<sup>k</sup>* corresponds to the total number of customers who rent at station *Si* during time *t*(*k*). The returning intensity of Poisson process *Returni*,*<sup>k</sup>* corresponds to the total number of returning customers. *Renti*,*<sup>k</sup>* and *Returni*,*<sup>k</sup>* can be calculated as follows:

$$Rent\_{i,k} = \lambda\_i(k) \times T \tag{8}$$

$$Return\_{i,k} = \mu\_i(k) \times T \tag{9}$$

where *λi*(*k*) and *μi*(*k*) correspond to the customer arrival rate for renting and returning respectively. *T* is the length of time slice *k*. It is worth noting that we divide a day into 24 time slices. For instance, *λi*(8) means the arrival rate for renting from 7:00 am to 8:00 am.

According to the *Renti*,*<sup>k</sup>* and *Returni*,*k*, the probabilities of renting (*Prent*) and returning (*Preturn*) at each step at time slice *k* are defined respectively as follows:

$$P\_{nnt} = \frac{Rent\_{i,k}}{Rent\_{i,k} + Renturn\_{i,k}} \tag{10}$$

$$P\_{return} = \frac{Return\_{i,k}}{Return\_{i,k} + Return\_{i,k}} \tag{11}$$

In Equations (10) and (11), we ignore the influence of customers' arrival order for renting and returning and only focus on the total number of customers that have arrived in this time period [20].

The above discussions indicate that the proposed RWTAW algorithm can describe the changing state of a station during a certain period. Algorithm 1 gives the details of calculating *OPi*,*k*. This proposed algorithm can be computed in O(*n*) operations. The total runtime of the algorithm is determined by the capacity of the station. Next, we will introduce how to choose the optimal station at which the system incentivizes the customer to rent bikes.



#### *4.2. The Optimal Station Selecting*

When a customer chooses his target station *Si* at time *t* on the APP (application program) of BSS, the system would incentivize him to rent bikes at another optimal station *OSt* if *Si* is not in its optimal state (*mi*(*t*) < *OPi*,*k*). Two issues need to be solved before deciding the optimal station: (1) The selected station should be as close as possible to the customer's target station, because customer may be unwilling to walk too far; and (2) The incentive measure should have a minimal impact on the station's sustainable working or can help the station be closer to its optimal state. Aiming to solve these two issues, we proposed the TLTF algorithm to determine the optimal station *OSt* which satisfies:

$$\text{OS}\_{l} = \underset{\text{S}\_{j}}{\text{arg}\,\text{max}} \,\delta\_{j,k} \,\forall \text{S}\_{j} \in \mathbb{C}\_{i} \tag{12}$$

where *Ci* is the set of candidate stations. We choose *Ci* based on the maximum distance *r* that customers are willing to walk [11]. Specifically, we select the stations which are no more than *r* from *Si* as the candidate stations in *Ci*:

$$\mathcal{C}\_{i} = \{ \mathcal{S}\_{\bar{\jmath}} \mid d(\mathcal{S}\_{i}, \mathcal{S}\_{\bar{\jmath}}) \le r \} \tag{13}$$

In Equation (12), *δi*,*<sup>k</sup>* represents the difference between the current state of station *Si* and its optimal state at time slice *k*, which is the key influence on *OSt*. Its definition is given as follows:

$$
\delta\_{i,k} = m\_i(t) - OP\_{i,k} \tag{14}
$$

If *δi*,*<sup>k</sup>* > 0, the rebalancing strategy will not be triggered when a customer chooses *Si* for renting. If *δi*,*<sup>k</sup>* -0, the system will incentivize him to rent bikes at station *OSt*. As defined in Equation (12),

the station *OSt* is the one which has the largest difference between its current state and the optimal state. The procedure to determine the optimal station *OSt* is described in Algorithm 2. The complexity of the algorithm is O(*n*).

#### **Algorithm 2:** TLTF


#### **5. Experimental**

#### *5.1. Data Set Description*

The datasets used in this paper were obtained from Bay Area, a bike sharing system in San Francisco. We have uploaded the datasets to github (https://github.com/TwinkleBill/babs\_open\_ data\_year\_3). The datasets include 68 stations' records from 1 September 2015 to 31 August 2016. Each record includes the status of station and customer's trajectory information. We give a visualization of the stations in Figure 3 and show the statistical information of the datasets in Table 2.

**Figure 3.** Map of the Bay Area bike sharing system. (**a**) Stations in Bay Area; (**b**) Stations in our simulation.

**Table 2.** Details of the Dataset.


#### *5.2. Performance Metrics for Bike-Sharing System*

Here we define three metrics to evaluate the performance of the rebalancing strategies in a bike sharing system.


$$\gamma\_{\text{lost}} = \frac{n\_{\text{lost}}}{n\_{\text{lost}} + n\_{\text{ridc}}} \tag{15}$$

where *nlost* is the number of customers who fail to rent or return a bike, *nride* is the number of customers who successfully ride a bike. The smaller the value of *γlost*, the better the performance of the system.

• The times of reposition: The number of bikes that need to be rebalanced (the truck-based rebalancing strategies) or the number of customers (the customer-oriented rebalancing strategies). In terms of the truck-based rebalancing strategies, the times of relocating bikes corresponds to the truck operating costs. For customer-oriented rebalancing strategies, the incentive cost also increases with the increasing number of incentive customers. To simplify things, the cost of delivering a bike is seen as the same as incentivizing a customer. Obviously, the smaller the times of reposition, the less the costs of rebalancing.

#### *5.3. Parameter Setting*

In the simulative bike sharing system, we need to set the following parameters: (1) *λ<sup>i</sup>* and *μi*, the customer arrival rates for renting and returning at station *Si* respectively. (2) The distribution parameters of trip duration *Ti*,*j*. (3) The walking distance *r* in the TLTF algorithm.

As discussed before, for each station in a BSS, the arrival processes of customers for renting and returning can be described by the Poisson process [43], which indicates that the time interval between two adjacent events should obey the exponential distribution. We analyzed the time interval distribution of the experimental data and found that most of the stations followed the exponential distribution. We take the historical records of the station named "San Jose Diridon Caltrain Station" from 7:00 a.m. to 8:00 a.m. (the morning peak [50], as shown in Figure 4) as an example. Figure 5a shows the time interval distribution of the renting process. It follows the exponential distribution with the parameter *<sup>λ</sup>* <sup>=</sup> <sup>1</sup> 11.5208. Thus, we utilize the Poisson process with *<sup>λ</sup><sup>i</sup>* <sup>=</sup> <sup>1</sup> 11.5208 to simulate the renting process of customers at this station from 7:00 a.m. to 8:00 a.m. Similarly, as shown in Figure 5b, the Poisson process with *<sup>μ</sup><sup>i</sup>* <sup>=</sup> <sup>1</sup> 8.4724 is used to simulate the returning process. In the RWTAW algorithm, the values of *λ<sup>i</sup>* and *μ<sup>i</sup>* of each station are different.

**Figure 4.** The renting and returning bikes by hours. (**a**) Renting by hours; (**b**) Returning by hours.

**Figure 5.** Statistical information of time interval and trip duration. (**a**) Rent; (**b**) Return; (**c**) Trip Duration.

The trip duration is also an important parameter in the simulation system. As shown in Figure 5c, the distribution of trip duration obeys the lognormal distribution (the mean is 6.2806 and the sigma is 0.7032). In our simulation system, the trip duration of each trip is generated by this lognormal distribution.

In order to demonstrate the effectiveness of our simulated bike sharing system, the result of one simulation of the station named "San Jose Diridon Caltrain Station" is illustrated in Figure 6. We plot the variation of the number of bikes rented by customers. The changing pattern of the simulation is in line with that of the empirical data. According to the results of the simulation, we have confidence that the simulated bike sharing system with the aforementioned parameters can reproduce the results of a real system.

In TLTF, the value of walking distance *r* is crucial, because it will take more incentivizing bonus to encourage customers to rent bikes at stations further away. Furthermore, the *r* is also related to the computational cost of the algorithm. In other words, a large value of *r* means lots of stations will be selected as candidate stations by the algorithm when it calculates the *OSt*. On the other side, the increase of *r* may help to select the optimal station for renting. Figure 7 shows the influence of distance *r* on the performance of the system. When *r* increases from 600 m to 1000 m, the unworking time and the proportion of lost customers decrease. Because the candidate stations will include more stations with an increase of the walking distance *r*. For example, in the datasets, *S*<sup>14</sup> is the only candidate of *S*<sup>2</sup> when *r* = 600 m, while the candidate stations consist of *S*4, *S*<sup>5</sup> and *S*<sup>14</sup> when *r* = 1000 m. However, the result will be bad if the walking distance is too far. As we can see, the unworking time and the proportion of lost customers increase when *r* = 1200 m. Although a large *r* increases the number of candidate stations, it also leads to an increase in the size of the intersection of candidate stations. These stations will be recommended by the system. In this case, they will become empty

or full soon. So the unworking time and the proportion of lost customers increase when *r* = 1200 m. We also care about the relationship between the probability of a customer's acceptance of walking distance and people's travel behavior [51], and count each station's walking distance to its adjacent stations which are not more than 1000 m away. We find that the stations uniformly distribute within 1000 m away from a center station. We think the acceptance probability of this walking distance is also suitable for customers, because they have a wide choice for purpose station [11]. Therefore, we set *r* = 1000 m.

**Figure 6.** Simulation Results of station named "San Jose Diridon Caltrain Station".

**Figure 7.** System performance with different *r* in TLTF.

In addition, we also analyze bike usage on workdays and weekends. As shown in Figure 8, the patterns of the system on workdays are similar, that is, two peaks respectively appear at 10:00 a.m. and 18:00 p.m. One is the morning peak, the other is the evening peak (the two peaks for whole system are 8:00 a.m. and 18:00 p.m., as shown in Figure 4a). However, the patterns on weekends are distinctly different. On weekends, the usage is lower than that on workdays, and is irregular. Such a phenomenon is consistent with what we have mentioned before, that is, commuting time is an important factor which leads to the stable traffic pattern on workdays. Therefore, in this work, we only simulate the operation of the bike sharing system on workdays. The parameters of the system are the same on each workday.

Note that in the dataset of Bay Area, stations are distributed into four districts (as shown in Figure 3a). There are no trajectories between any two different districts. The reason is that customers generally do

not choose bikes for long-distance travel. Therefore, stations in our simulation system are all in "San Jose" which has 15 stations (as shown in Figure 3b). The probability of choosing the returning station for the customer satisfies the uniform distribution. Experiments indicate that such an assumption does not influence the time interval distribution of a customer's arrival process for returning to a station, which still fits well with the exponential distribution.

**Figure 8.** Bike usage on workdays and weekends. (**a**) Monday; (**b**) Tuesday; (**c**) Wednesday; (**d**) Thursday; (**e**) Friday; (**f**) Saturday; (**g**) Sunday.

#### *5.4. Results*

In order to validate the proposed strategy, we compared it with three scenarios: the system with no rebalancing, the static rebalancing strategy and the dynamic rebalancing strategy. Note that our target is to verify the effectiveness of the rebalancing strategy, so we do not take too much care on the truck routing problem and the price strategy of incentive in this paper.


Figure 9 shows the system's performance with different rebalancing strategies. The system with no rebalancing gets the worst performance on the unworking time and the proportion of lost customers. Compared with the NR strategy, the SR strategy decreases by nearly 40% and 39% on these two metrics, respectively. With the decrease of the scheduling time interval, the effectiveness of the dynamic strategy improves (in *DRi*, *i* = 1, 4, 8). As shown in Figure 9, among the three DR strategies, the one which has the shortest scheduling time interval (*i* = 1) has the least unworking time and the smallest proportion of lost customers. However, its times of reposition is the largest. According to the simulation results, the proposed rebalancing strategy can obtain the smallest proportion of lost customers with the proper times of reposition. At the same time, the unworking time is also within our acceptance.

**Figure 9.** Comparison of system performance.

In BSS, the traffic pattern changes rapidly. The rapid response of the system is important for offering customers a better service. Thus, the rebalancing algorithm should be in real-time. For static strategies, they are implemented at the system's idle time, which means there is no user in the system during this period. So we do not have to pay attention to the real-time performance of static strategies. Table 3 shows the runtime of our strategy and the DR strategy [20] to calculate a station's optimal state at a specific hour. There are four stations. The capacity of each station is 11, 15, 19 and 27 respectively. As we can see, the increase of the station's capacity increases the computing time. The time cost of the DR strategy is greater than ours. The results show that our strategy is real-time and efficient.


**Table 3.** Comparison of runtime of calculating the optimal state.

#### **6. Conclusions**

In this paper, we designed a customer-oriented rebalancing strategy for bike sharing systems. On the basis of the one-dimension Random Walk process with two absorption walls, we infer the probability that a station becomes empty or full, which is useful to both system operators and customers. The proposed strategy includes two algorithms called RWTAW and TLTF respectively. The RWTAW calculates the optimal state of the station. The TLTF determines the optimal station at which the system encourages customers to rent. Compared with the truck-based static and dynamic rebalancing methods, our strategy has the best performance on decreasing the imbalance of the system while keeping the rebalancing cost as low as possible. In future work, we will explore how to predict the number of bikes for each station to improve customer satisfaction.

**Author Contributions:** conceptualization, P.Y. and F.H.; methodology, P.Y.; software, P.Y.; validation, P.Y., F.H. and J.P.; formal analysis, P.Y.; investigation, F.H.; resources, J.P.; data curation, P.Y.; writing–original draft preparation, P.Y. and F.H.; writing–review and editing, P.Y., F.H. and J.P.; visualization, P.Y.; supervision, P.Y.; project administration, J.P.; funding acquisition, J.P.

**Funding:** This research was funded by the National Key Research and Development Project with Grant No. 2017YFB0202403, Sichuan Key Research and Development Program (2017GZDZX0003-02, 2019KJT0015-2018GZ0098), and the Sichuan University-Zigong Science and Technology Fund (2018CDZG-15).

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

#### **References**


c 2019 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/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Energies* Editorial Office E-mail: energies@mdpi.com www.mdpi.com/journal/energies

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18