DE GRUYTER OPEN HUNG. J. IND. CHEM. 2017 45(1), 1–84 Table of Contents Dynamical Temperature from the Phase Space Trajectory ANDRÁS BARANYAI ............................................................................................................................................ 1-3 Separation of Gases by Membranes: The Effects of Pollutants on the Stability of Membranes NÁNDOR NEMESTÓTHY ...................................................................................................................................... 5-8 Selective Removal of Hydrogen Sulphide from Industrial Gas Mixtures Using Zeolite NaA TAMÁS KRISTÓF .............................................................................................................................................. 9-15 Dynamic Simulator-Based Apc Design For A Naphtha Redistillation Column LÁSZLÓ SZABÓ, KLÁRA KUBOVICSNÉ STOCZ, LAURA SZABÓ, SÁNDOR NÉMETH, AND FERENC SZEIFERT .......... 17-22 Replacement of Biased Estimators with Unbiased Ones in the Case of Student's t-Distribution and Geary’s Kurtosis GERGELY TÓTH AND PÁL SZEPESVÁRY ........................................................................................................... 23-27 Formation, Photophysics, Photochemistry and Quantum Chemistry of the Out-Of-Plane Metalloporphyrins ZSOLT VALICSEK, MELITTA P. KISS, MELINDA A. FODOR, MUHAMMAD IMRAN, AND OTTÓ HORVÁTH .................. 29-36 Improvement Possibilities of Heterogeneous Photocatalysis with the Aim of In-Field Use ORSOLYA FÓNAGY, PÉTER HEGEDŰS, ERZSÉBET SZABÓ-BÁRDOS, ANNAMÁRIA DOBRÁDI, AND OTTÓ HORVÁTH ...................................................................................................................................................... 37-48 Adapting the SDEWES Index to Two Hungarian Cities VIKTOR SEBESTYÉN, VIOLA SOMOGYI, SZANDRA SZŐKE, AND ANETT UTASI ..................................................... 49-59 Surface Energy Heterogeneity Profiles of Carbon Nanotubes with a Copolymer-Modified Surface Using Surface Energy Mapping by Inverse Gas Chromatography FRUZSINA GERENCSÉR, NORBERT RIEDER, CSILLA VARGA, JENŐ HANCSÓK, AND ANDRÁS DALLOS .................. 61-66 Initial Electrical Parameter Validation in Lead-Acid Battery Model Used for State Estimation BENCE CSOMÓS, DÉNES FODOR, AND GÁBOR KOHLRUSZ ............................................................................... 67-71 Simulating Ion Transport with the NP+LEMC Method. Applications to Ion Channels and Nanopores. DÁVID FERTIG, ESZTER MÁDAI, MÓNIKA VALISKÓ, AND DEZSŐ BODA ............................................................... 73-84 . László Rácz HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) Editorial preface hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0021 IN MEMORIAM JÁNOS LISZI This issue is dedicated to Professor János Liszi who was the Rector of the University of Veszprém in 1989-1995. The University acquired this name during his term in an era of severe political and economic challenges (today: University of Pannonia). Under his leadership the University took on these challenges and has been modernized to become an institution with a wider scope. He was the head of the Department of Physical Chemistry for many years. The Department was also reshaped during this period. His early death was unexpected and mourned. He influenced the life of many and left his mark on his Alma Mater in Veszprém, but he also had an impact on other institutions as the commemoration of László Rácz below indicates. László Rácz was the Head of the Chemistry Department in Eger when Professor Liszi joined to teach Physical Chemistry there. It was an honor to organize and edit this issue to salute my respected supervisor. Veszprém, 10 October 2017. Dezső Boda Issue Editor THE EGER-CONNECTION OF PROFESSOR LISZI He was a tough and righteous man towards himself just as towards his colleagues in Veszprém. I can say that the same was true regarding his relationship with the colleagues and students of the Chemistry Department at the Eszterházy Károly College in Eger. That is why we got to like him so much. When a position for teaching Physical Chemistry became available in our Department, it was not easy to find the right person to fill the boots left by an excellent colleague who had held this position for 30 years. According to the urgent need to yield the proper knowledge to our students, our director indicated to Professor Liszi that there was a problem to be solved in Eger. He did not hesitate and drove 230 km to Eger right away. He assisted with our teaching work for many years. A fruitful professional relationship had been established between Veszprém and Eger. We gladly accepted the support offered by Professor Liszi. Beyond his Physical Chemistry courses in the curriculum, he regularly gave seminars about various topics not only to our students but also to colleagues and students from other departments. An example that I witnessed characterizes his professional work. After a class, a student from the Physics programme approached him and said, “I will enter the Chemistry programme because this was so beautiful and logical”. And he did. He became a Physics-Chemistry teacher and now works in the Physics Department. Since 2015, our College has functioned as a University. I am sure that the support of Professor Liszi has contributed to the development of our institution considerably. He has also contributed to my personal development because he encouraged me to defend my habilitation thesis. He also contributed to the education of our team beyond the field of Physical Chemistry. During excursions to Tihany and at The Benedictine Order he acted as a guide and demonstrated his vast knowledge of History and classical culture. The expertise of our Department, making wine, also belongs to classical culture. It was a pleasure to host him and his wife for wine tastings on several occasions. At such times, we discussed the science of wine-making. Time flew and we got onto the subject of the thermodynamics of multicomponent homogeneous systems after a few glasses of wine. Sadly the Blessing of St. John had to come early because he had an early start the following morning. The minimum 1000 m session in the swimming pool could not be left unperformed. Together with him we profess that technical details can only come after the firm scientific basics have been laid down. He achieved a lot for higher education and research in Eger. We are grateful for his efforts We respected and loved him not only for his knowledge, but also for his humanity. Eger, 12 October 2017 HUNGARIAN JOURNAL
OF INDUSTRY AND CHEMISTRY Vol. pp.45(1) pp. 1–3 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0001 DYNAMICAL TEMPERATURE FROM THE PHASE SPACE TRAJECTORY ANDRÁS BARANYAI * Department of Theoretical Chemistry, Eötvös Loránd University, P. O. Box 32, Budapest, 1518, HUNGARY In classical statistical mechanics the trajectory in phase space represents the propagation of a classical Hamil- tonian system. While trajectories play a key role in chaotic system theory, exploitation of a single trajectory has yet to be considered. This work shows that for ergodic dynamical systems the dynamical temperature can be derived using phase space trajectories. Keywords: temperature, statistical mechanics A Hamiltonian system of N mass points is related to the microcanonical ensemble of statistical mechanics when the volume, V, and the energy, E, of the system are fixed. A microscopic state of the system is a point in phase space represented by a 6N-dimensional vector, ),...,,,...,( 2121 NN pppqqqΓ  . The time evolution of this point is described as a trajectory in phase space. In microcanonical ensembles, since the energy is fixed, the evolution of the system is restricted to a 6N-1 dimensional hypersurface. The logarithm of the area of this hypersurface yields the entropy, S, of the system. Since the entropy of a stable thermodynamic system is a monotonically increasing function of the internal energy, the accessible phase space increases with energy. This provides the opportunity to determine the thermodynamic derivative, TES NV /1)/( ,  , where T is the absolute temperature. The trajectories within the limit of infinite time perfectly cover the 6N-1 dimensional hypersurface according to the ergodic hypothesis. This means that the length of the trajectory for a period of time,  , is proportional to the microcanonical entropy: .limln 0 2             ΓdtkS (1) Eq.(1) is a line integral measuring the length of the path of the moving system in phase space. This equation can be rewritten as: .limlnlimln 00 2                                   ΓΓ H dtk H dtkS (2) *Correspondence: baranyai@chem.elte.hu Let us consider two Hamiltonian systems, H and H’, with only a minor difference between their energy parameters, E and E’. The distance covered in phase space by these systems is different because the two energies belong to different areas of the hypersurface to be covered by the trajectories. If the difference in entropies between the two systems is calculated, the quotient in the argument of the logarithm makes it possible to use equality instead of a mere proportio- nality. . ' ln ' ln lim lim ln)()( 0 0                                                                             Τ Γ Τ Γ Γ Γ H H k H H k H dt H dt kESES       (3) In Eq.(3) the time integral was replaced with the product of time and the time-average of the integrands. The Hamiltonian, H’, can be approximated as .Γ Γ d H HH    (4) Then . 2 2 Γ ΓΓΓ d HHH         (5) BARANYAI Hungarian Journal of Industry and Chemistry 2 Using Eq.(5), Eq.(3) can be rewritten as . ln)()( 2 2 2 2 Γ Γ Γ Γ Γ ΓΓ                             H d H k H d HH kESES (6) In the second equality of Eq.(6), the xx  )1ln( app- roximation was exploited. To obtain the temperature, division with the ener- gy difference should be included: . 1 )()( 2 2 2 2 2 TH H k d H H d H k EE ESES                 Γ Γ Γ Γ Γ Γ Γ (7) The dynamical temperature was derived first by Rugh as a time average of )/( 2 HH  on the energy surface [1]. . 1 2 H H kT    (8) Later, Butler et al. applied a similar approach and derived Eq.(8) only for the configuration part of the phase space and also checked the performance of the method numerically [2]. The essence of both derivations was to transform the phase-space vector, Γ , from system E into Γ of system E’ by a vector containing the Hamiltonian gradient. In this way arelationship was established between the phase-space points of the two systems. Our result uses the average speed of evolution of trajectories in phase space, therefore, a Hamiltonian gradient is not required to connect the two sets of phase- space points. Expansion of the energy in Eq.(4) is sufficient to relate the two trajectories. Certainly, the results are identical which proves the validity of this alternative approach while, at the same time, it is a nice example of ergodicity: the method which uses the difference in area of hypersurfaces in an ensemble- average fashion leading to results that are identical to the method that uses relative trajectory propagation in time. The explicit form of Eq.(7) can be written as , / /3 1 1 22 2 1        N i i N i m m kT q q p (9) where  is the position-dependent potential in the Hamiltonian and m is the mass of a particle. Eq.(9) contains two terms that feature both in the numerator and denominator. The quotient of the first terms is the well-known kinetic temperature and the quotient of the second terms is the configurational temperature. In actual calculations these two temperatures can be calculated contrary to the full form of Eq.(9) which contains dimensional discrepancies in both the numerator and denominator. The configurational temperature calculation was suggested as an algorithmic check for Monte Carlo computer simulations where the Boltzmann temperature is an input parameter [2]. For completeness, the third, very simple method of deriving the configurational temperature is mentioned [3].The previous derivations did not state whether the two temperatures, the kinetic and configurational, are equal. Using a trivial derivation it is shown that the two quotients are equal [3]. A further advantage of this method is that there is no need for the condition of N to accept the validity of a microcanonical result in terms of the canonical ensemble. The temperature from the momenta of the particles is . 3 1 1 2   N i i Nkm T p (10) If this temperature remains constant over time: ,0 11    N i ii N i ii CCT Fppp  (11) where, for the sake of simplicity, the constant factor is denoted in front of the sum as C and  ii qF  / is the Newtonian force, it is also possible to write that .0 )( 1 ,, 2 2 1                    N i zyx ii i N i iiii m p FC CT     F FpFp  (12) At equilibrium in isotropic systems the different Cartesian directions are equivalent and there is no correlation between velocities and position-dependent quantities. Thus, DYNAMICAL TEMPERATURE FROM THE PHASE SPACE TRAJECTORY 45(1) pp. 1–3 (2017) 3 . 3 2 2 2    q q kT p m (13) Eq.(13) can be written in this form when T is mani- pulated. The results would be the same. Therefore, if the expectation value of the temperature is constant, the kinetic and configurational temperatures are equivalent. As for the temperature derivation from the phase- space trajectory it is important to note that the  condition is essential. Finite segments of trajectories do not cover the hypersurface completely and only contain minimal information about the temperature and force distributions of the system. REFERENCES [1] Rugh, H.H.: Dynamical approach to temperature, Phys. Rev. Lett., 1997 78(5), 772-774 DOI:10.1103/PhysRevLett.78.772 [2] Butler, B.D.;Ayton, G.A.;Jepps, O.G.;Evans, D.J.: Configurational temperature: Verification of Monte Carlo simulations, J. Chem. Phys., 1998 109, 6519- 6522 DOI:10.1063/1.477301 [3] Baranyai, A.: On the configurational temperature of simple fluids, J. Chem. Phys., 2000 112, 3964- 3966 DOI:10.1063/1.480995 HUNGARIAN JOURNAL
OF INDUSTRY AND CHEMISTRY Vol. pp.45(1) pp. 5–8 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0002 SEPARATION OF GASES BY MEMBRANES: THE EFFECTS OF POLLUTANTS ON THE STABILITY OF MEMBRANES NÁNDOR NEMESTÓTHY* Research Institute of Bioengineering, Membrane Technology and Energetics, University of Pannonia, Egyetem u. 10, Veszprém, H-8200, HUNGARY The long-term stability of membranes is determined mainly by their sensitivity to pollutants. Their stability was tested using a novel, multichannel measuring system, which is based on pressure differences. This measuring system is suitable to determine the changes in permeability of polymer membranes. The damaging effects of H2S, BTX and n-dodecane were investigated in terms of polyimide gas separation membranes using nitrogen gas. Keywords: multichannel test equipment, pressure difference, hydrogen sulphide, BTX 1. Introduction Previously, it was thought that the stability of memb- ranes is determined by the mechanical stress (shear) and the natural aging of polymers. Recently, however, it has been confirmed that their stability is limited mainly by the sensitivity of membranes to certain pollutants. These aggressive compounds, pollutants, e.g. chlorine, hydrogen sulphide, hydrocarbons, etc., may damage the structure of the polymer, thus its physical and chemical properties change, and consequently the permselectivity of the membrane changes, as well. It is known that polymeric reverse osmosis memb- ranes are sensitive to strong oxidising agents, especially chlorine compounds [1], therefore, intensive research has been conducted to avoid or at least reduce any damage [2-3]. Similar levels of membrane degradation are observed in proton-exchange membrane (PEM) fuel cells and batteries containing membranes, where oxi- dising compounds are in contact with the membranes, as well [4-5]. The stability of polymeric gas separation memb- ranes has hardly been investigated. The long-term eff- ects of H2S on inorganic membranes has been studied by Australian researchers at low concentrations (50 ppm) [6]. However, H2S may not only cause long-term, but immediate damage, mainly in the form of swelling, which strongly influences the gas transport properties of membranes. Koros and co-workers presented the effects of extremely high H2S concentrations on a polymer membrane (50,000 - 100,000 ppm) [7]. In the field of membrane technology, sensitivity can be measured by a sort of effectiveness unit. The *Correspondence: nemesn@almos.uni-pannon.hu product of the concentration and time period yields a value where the effectiveness of the membrane decr- eases from 95 to 90 and then to 70% of its original value (e.g. 1000 ppm*hour means that the membrane was exposed to 1000 ppm of pollutant for 1 hour, or 0.1 ppm of pollutant for 10,000 hours during the tests). According to the literature these types of measurements have yet to be published for gas separation membranes, thus the aim of this work was to design, construct and operate a piece of test equipment that conducts reliable laboratory tests. For the determination of stability, direct and indirect methods can be used to measure the gas volumes passing through the membranes. Direct methods are usually preferred, and – if the composition of the gas is known – are more exact than indirect ones. However, when the gas composition varies and small amounts of gases need to be measured, indirect methods are often more suitable. In this work an indirect method based on a pressure differential technique was chosen, where the pressure of a closed vessel is measured and the varying pressure yields information about the volume of the gas passing through the membrane. During the investigation the effects of pollutants on the permeability of nitrogen was to be studied. The following pollutants were used:  compounds containing sulphur at associated gas  BTX mixture (benzene, toluene, xylene)  heavy hydrocarbons In this research the aim was to determine quantitatively the effects of pollutants on the membranes to define the tolerance range of particular membranes. NEMESTÓTHY Hungarian Journal of Industry and Chemistry 6 Figure 1. The small modules constructed for the tests 2. Experimental For this series of measurements, polyimide gas separation membranes (synthesised by UBE) were used. They were taken from a hollow fibre module and can accurately model the properties of industrial gas separation membranes. From the hollow fibres small modules containing 6 capillaries were constructed (Fig.1) and their ends were closed, thus their tests were carried out in a “sack” configuration. In the design of the test system it was important that several parallel measurements should be conducted and the measuring channels combined with each other. The scheme of a measuring channel can be seen in Fig.2 The gas was introduced into the measuring system through valve V1 (which can be adjusted by valve V2 if necessary). Before measurements were taken the pressure of the vessel was checked by pressure transducer PT1. To start the test the pressure was adjusted by regulating valve PV1, which was checked by pressure transducer PT2. Then the membrane was installed into the thermostatic system in a way that ensured its mobility was not restricted, thus the permeation of gas could not influence the flux. A photograph of the measuring system is shown in Fig.3. For the permeability measurements nitrogen gas from a cylinder was used (99.5%; Messer Hungarogáz Kft., Hungary). The permeability of the membrane was determined from the pressure of the vessel and the transmembrane pressure measured on-line during the experiments. In the experiments, H2S, methyl mercaptan and ammonia (compounds containing S or N), oleic acid, ethyl Figure 2. The scheme of the test system Figure 3. The test system alcohol, moreover, a benzene-toluene-xylene mixture (BTX) and n-dodecane (as a heavier hydrocarbon) were used as pollutants. For the stability experiments the small membrane modules were put in a closed vessel (Fig.4) where the headspace was saturated with the given pollutant. The vessels were placed in a thermostatic incubator at 27 °C usually for between 1 and 7 days. Certain materials (e.g. BTX) damaged the epoxy resin glue used to adhere the fibres of the membranes, therefore, these experiments were repeated using poly- ether-sulfone glue instead. 3. Results The nitrogen permeability of the membranes was determined before and after the incubations. In the preliminary experiments, the ammonia solution and methyl mercaptan severely damaged the surface of the membranes, thus no flux could be measured. Oleic acid and ethyl alcohol hardly influenced the flux, while BTX, H2S and n-dodecane changed the permeability of nitrogen considerably. For further investigation of the pollutants an experimental design was constructed, using appropriate Figure 4: The membranes in the closed vessel THE EFFECTS OF POLLUTANTS ON THE STABILITY OF MEMBRANES 45(1) pp. 5–8 (2017) 7 Table 1. The parameters of the experimental design pollu- tant cmin ppm cav ppm cmax ppm tmin d tav d tmax d H2S 100,000 300,000 500,000 1 3.5 7 BTX mixt. 1,000 750 500 1 3.5 7 dode- cane 1,000 5,500 10,000 1 3.5 7 statistical methods. The parameters selected were the concentrations of the pollutants (minimum, maximum and average) and the incubation time (minimum, maximum and average). The Statistica 8 computer program was applied to the design presented in Table 1. Firstly the effect of H2S on the nitrogen permeability of the membranes was measured. The experimental results are presented in Fig.5. It can be seen that the permeability of nitrogen increased even when the concentrations of pollutants were low and rose by using higher concentrations and longer periods of exposure. From the H2S concentration and the incubation time it was possible to calculate a special parameter of exposure with the unit of ppm*h. Permeability was presented as a function of this parameter (Fig.6), where an almost linear relationship was observed. The results suggest that the process can be described as a first-order reaction, which means that no safe limit can be determined where H2S is regarded as harmless, on the contrary, it should be considered at all times. The effect of the BTX mixture was studied using a similar methodology. The results are summarized in Fig.7. This figure shows that the permeability of nitrogen increased at low concentrations of the BTX mixture. At higher concentrations and over longer periods of time, no further significant changes were observed. The exposure parameter was also calculated in the unit of ppm*h and permeability was presented as its function (Fig.8). The diagram can be described as a saturation-type curve. 2**(2-0) design; MS Residual=49886.27 > 2000 < 2000 < 1500 < 1000 < 500 < 0 Figure 5. The effect of H2S on the nitrogen permeability Exposure [M ppm * h] 0 20 40 60 80 100 P e rm e a b ili ty c h a n g e [ % ] 0 500 1000 1500 2000 2500 Figure 6. Permeability changes against exposure to H2S 2**(2-0) design; MS Residual=2675.592 > 400 < 360 < 310 < 260 < 210 < 160 < 110 < 60 Figure 7: The effect of the BTX mixture In the last series of experiments, the effect of exposure to n-dodecane was investigated experi- mentally. The results are presented in Fig.9. n-dodecane caused – unlike H2S and BTX – a reduction in the permeability of nitrogen even at low concentrations. The flux fell to zero at higher concentrations and over longer incubation periods. Permeability was investigated as a function of exposure (Fig.10). The process can be described as a first-order reaction, thus the effect of heavier hydrocarbons, e.g. n-dodecane, should always be considered. Exposure [k ppm * h] 0 20 40 60 80 100 120 140 160 180 p e rm e a b ili ty c h a n g e [ % ] 0 100 200 300 400 Figure 8: Permeability changes against BTX exposure NEMESTÓTHY Hungarian Journal of Industry and Chemistry 8 2**(2-0) design; MS Residual=40.32141 > 100 < 96 < 76 < 56 < 36 Figure 9: The effect of n-dodecane 4. Conclusion The long-term stability of polyimide gas separation membranes was tested against various pollutants: H2S, a BTX mixture and n-dodecane. These compounds significantly affected the nitrogen permeability of the membranes which were described by using a special parameter of exposure. It was found that H2S and the BTX mixture increased the permeability, while n- dodecane reduced the permeability of the membranes. Further investigations are planned to investigate the effect of other pollutants, moreover, to determine the permeability of additional gases, e.g. carbon dioxide, methane, etc. Acknowledgement We acknowledge the financial support of Széchenyi 2020 under the EFOP-3.6.1-16-2016-00015 project. This research was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences. REFERENCES [1] Glater, J.; Hong, S.K.; Elimelech, M.: The search for a chlorine-resistant reverse osmosis membrane, Desalination, 1994 95(3), 325-345 DOI: 10.1016/0011-9164(94)00068-9 Exposure [ k ppm * h] 0 50 100 150 200 p e rm e a b ili ty c h a n g e [ % ] 0 10 20 30 40 50 60 70 Figure 10. Change in permeability against exposure to n-dodecane [2] Zhang, Y.; Zhao, C.; Yan, H.; Pan, G.; Guo, M.; Na, H.; Liu, Y.: Highly chlorine-resistant multilayer reverse osmosis membranes based on sulfonated poly(arylene ether sulfone) and poly(vinyl alcohol), Desalination, 2014 336, 58-63 DOI: 10.1016/j.desal.2013.12.034 [3] Gohil, J.M.; Suresh, A.K.: Chlorine attack on reverse osmosis membranes: Mechanisms and mitigation strategies, J. Membr. Sci., 2017 541, 108-126 DOI: 10.1016/j.memsci.2017.06.092 [4] Huang, X.; Pu, Y.; Zhou, Y.; Zhang, Y.; Zhang, H.: In-situ and ex-situ degradation of sulfonated polyimide membrane for vanadium redox flow battery application, J. Membr. Sci., 2017 526, 281- 292 DOI: 10.1016/j.memsci.2016.09.053 [5] Lapicque, F.; Belhadj, M.; Bonnet, C.; Pauchet, J.; Thomas, Y.: A critical review on gas diffusion micro and macroporous layers degradations for improved membrane fuel cell durability, J. Power Sources, 2016 336, 40-53 DOI: 10.1016/j.jpowsour.2016.10.037 [6] Uhlmann, D.; Smart, S.; Diniz da Costa, J.C.: H2S stability and separation performance of cobalt oxide silica membranes, J. Membr. Sci., 2011 380(1), 48-54 DOI: 1016/j.memsci.2011.06.025 [7] Kraftschik, B.; Koros, W.J.; Johnson, J.R.; Karvan, O.: Dense film polyimide membranes for aggressive sour gas feed separations, J. Membr. Sci., 2013 428, 608-619 DOI: 10.1016/j.memsci.2012.10.025 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 9–15 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0003 SELECTIVE REMOVAL OF HYDROGEN SULPHIDE FROM INDUSTRIAL GAS MIXTURES USING ZEOLITE NaA TAMÁS KRISTÓF* Department of Physical Chemistry, Institute of Chemistry, University of Pannonia, Egyetem u. 10., Veszprém, H-8200, HUNGARY Hydrogen sulphide removal from simple gas mixtures using a highly polar zeolite was studied by molecular sim- ulation. The equilibrium adsorption properties of hydrogen sulphide, hydrogen, methane and their mixtures on dehydrated zeolite NaA were computed by Grand Canonical Monte Carlo simulations. Existing all-atom intermo- lecular potential models were optimized to reproduce the adsorption isotherms of the pure substances. The ad- sorption results of the mixture, also confirmed by IAST calculations, showed very high selectivities of hydrogen sulphide to the investigated non-polar gases, predicting an outstanding performance of zeolite NaA in techno- logical applications that target hydrogen sulphide capture. Keywords: Hydrogen sulphide, Zeolite, Selectivity, Gas mixture, Molecular simulation 1. Introduction Hydrogen sulphide is a highly toxic, acidic and corrosive substance. It is present naturally in landfills, natural and biogases, as well as in several synthesis gases. One of its main anthropogenic sources is the processing of crude oil in industrial refineries, where hydrodesulfurization (HDS) of a variety of streams (e.g. engine fuels) produces hydrogen sulphide-containing gas mixtures, which need to be purified. The economic removal of hydrogen sulphide is a long-standing task of the oil and gas industry. Adsorptive separation involves the use of solid substrates with a specific affinity with particular compounds of the mixture. Zeolites have been applied as catalysts in the petrochemical industry for a relatively long time and these materials can also be used for purification/separation purposes. Zeolites are crystalline aluminosilicates consisting of a three- dimensional framework of SiO4 and AlO4 tetrahedra of a highly regular porous structure [1]. The typical size of zeolitic micropores is similar to that of many small molecules. In contrast to various other adsorbents, zeolites generally endure high temperatures and pressures well, and can tolerate harsh chemical environments. The Si/Al ratio is a key factor in the application of zeolites. Zeolites with lower Si/Al ratios are more hydrophilic, whereas high-silica zeolites often possess fewer structural defects. These latter adsorbents are preferred in the separation of non-polar gases. Zeolite NaA is a synthetic microporous zeolite which accommodates extraframework Na+ ions. It exhibits an especially high affinity with small polar *Correspondence: kristoft@almos.uni-pannon.hu molecules such as water. The adsorption and separation properties of zeolite NaA have already been examined in several experimental [2-6] and theoretical/simulation [7-16] works. In our laboratories, the adsorption characteristics of zeolite NaA and its performance as a drying agent by classical atomistic simulations [10, 12- 14] were studied, and new intermolecular potential models for this zeolite [12-14] proposed. Our models were optimized for the study of the selective adsorption of water from its mixtures with less polar or non-polar molecules like simple alcohols, carbon monoxide, hydrogen and methane. In this paper, the selective removal of hydrogen sulphide by zeolite NaA is investigated. Molecular simulation predictions for mixture adsorption from two- and three-component gas mixtures containing hydrogen sulphide (H2S), hydrogen (H2) and methane (CH4) are presented. 2. Models and simulation details Zeolite NaA [17-18] is of LTA framework type, the structure of which belongs to the Fm-3c space group with a lattice parameter of 2.4555 nm. The three- dimensional cubic arrangement of its framework atoms is comprised of three kinds of rings with four (4R), six (6R) or eight (8R) oxygen atoms. The interconnection of 4R and 6R rings forms nearly spherical cages (sodalite cages) and these cages are linked by oxygen bridges, shaping straight channels of supercages with a maximum diameter of about 1.2 nm. The standard type of zeolite NaA has a Si/Al ratio of 1. In the present study, the unit-cell composition of the standard type of zeolite NaA was chosen: it consists KRISTÓF Hungarian Journal of Industry and Chemistry 10 of 576 framework atoms, namely 96 silicon, 96 aluminium and 384 oxygen atoms. The framework atoms were fixed at the atomic positions measured in X- ray diffraction experiments [17] and the 96 non- framework Na+ ions were allowed to move. According to the Löwenstein rule that prohibits AlOAl linkages, each AlO4 tetrahedron of this framework is connected to a SiO4 tetrahedron. Realistic and rigid all-atom intermolecular potential models, consisting of Lennard-Jones and Coulombic interaction sites, were used in the simulations. In these models, the interaction sites were fixed at their experimental atomic positions and assigned their Lennard-Jones energy (ε) and size (σ) parameters, as well as point charges (q). For zeolite NaA, the model that was developed earlier was modified slightly [14] by adding weak Lennard-Jones interaction sites with realistic size parameters [16, 19] to the (originally) pure Coulombic silicon and aluminium atoms, thus preserving the dominant role of oxygen atoms in the dispersion interactions of the framework. For hydrogen sulphide, a rigid four-site model proposed recently by Shah et al. [20] was adopted, in which the location of the charge parameter of the sulphur atom is offset on the HSH angle bisector towards the hydrogen atoms. An OPLS-AA model [21] was used for the methane molecules, and its Lennard-Jones H parameters were also applied to the H2 molecules, with partial charges on the atomic sites and on the molecular centre of mass [22]. Table 1 lists the potential parameters of the above models. Instead of the generally accepted Lorentz-Berthelot combining rule, the unlike Lennard-Jones interactions were computed by the combining rule proposed by Waldman and Hagler [23]. This combining rule links the behaviour of the unlike energy parameter εij to the relative sizes of atoms i and j and yields somewhat smaller values for the parameters εij and ij when ii≠jj. Song et al. [24] found that the experimental thermodynamic properties of pure methane can be reproduced more accurately using the Waldman-Hagler combining rule. Gas adsorption simulations were carried out by the standard grand canonical Monte Carlo methodology [25]. Total pressure p and partial pressures in the gas phase were specified by the chemical potentials of the components; in a diluted gas, these can be calculated from the ideal gas law [12]. The long-range Coulomb interactions were treated with the Wolf method [26-27] using a convergence parameter of  = 2/rc and cutoff radius of rc = L/2 (L is the length of the simulation box). The simulations involved an equilibration period of 5×107 Monte Carlo moves and an averaging period of 2×108 moves, consisting of 70% molecular insertion/deletion and 30% molecular translation steps. Since the random insertion of molecules is unable to take into account the inaccessibility of the sodalite cages by multiatomic molecules (the physical diffusion pathways to them), creation of H2S and CH4 molecules inside these cages was blocked artificially by placing repulsive dummy atoms at the centres of the cages. As H2 molecules are sufficiently small to pass through the windows of the sodalite cages, the insertion of H2 molecules into these cages was permitted. In either case, the transition of molecules into sodalite cages via translational trial moves was not artificially prevented. In addition to the adsorption loading, the isosteric heat of adsorption was calculated using the equation: TV a a Tp b b n U n H q ,,                       , (1) where Hb and Ua stand for the residual enthalpy and residual internal energy, respectively, n is the number of moles of the substance in the adsorbed (a) or bulk (b) phases. In the grand canonical ensemble, the second part of the equation can be determined from the particle number fluctuations of the simulation and the cross- correlation of potential energy and particle number fluctuations [28-29]. Assuming the ideal gas adsorbates, the first part of the equation is equal to RT, where R is the gas constant and T is the temperature. Predictions for mixture adsorption were also made using the ideal adsorbed solution theory (IAST) [30-31], which is an analogue of the ideal Raoult’s Law. Using IAST, mixture adsorption loadings at a given T can be obtained from single-component adsorption loadings by determining the bulk pressure of each component po at the same spreading pressure  of the adsorbed phase: )(/ o i b i a i ppyy  . (2) Table 1. Lennard-Jones energy (ε), size parameters (σ) and partial charges (q) for the models used in this work (d is the bond length, k is the Boltzmann constant). interaction site σ / nm (ε/k ) / K q / electron charge position in the structure/molecule Na+ (NaA) 0.250 100 0.60 random positions in supercages Si (NaA) 0.230 22.0 2.40 experimental atomic positions [17] Al (NaA) 0.240 16.5 1.80 experimental atomic positions [17] O (NaA) 0.330 190 -1.20 experimental atomic positions [17] S (H2S) 0.360 122 - dS-H = 0.134 nm H (H2S) 0.250 50.0 0.21 HSH angle: 92° XS (H2S)* - - -0.42 dS-X = 0.03 nm C (CH4) 0.350 33.21 -0.24 dC-H = 0.109 nm H (CH4) 0.250 15.1 0.06 HCH angle: 109.47° H (H2) 0.250 15.1 0.4829 dH-H = 0.0741 nm centre of mass (H2) - - -0.9658 geometric centre of the linear H2 molecule *off-atom site on the H–S–H angle bisector towards the hydrogen atoms SELECTIVE REMOVAL OF HYDROGEN SULPHIDE USING ZEOLITE 45(1) pp. 9–15 (2017) 11 Here, y i is the mole fraction of component i, and )(o ip is given implicitly by:  o 0 o ln)()( ip a ii pdpn A RT p , (3) where A is the surface area of the adsorbent. 3. Results and discussion The intermolecular potential models were tested by determining equilibrium adsorption isotherms for pure H2S, H2 and CH4 on zeolite NaA. Experimental data at 298 K are available for H2S [32] and H2 [8] and nearly room-temperature (T = 283 K) data for CH4 were taken from [33]. Fig.1 shows that the models are largely able to reproduce the experimental adsorption data for these substances. The reproduction of the experimental isotherm is quite good for H2S and H2. In the case of CH4, the extent of overestimation of the experimental results at 283 K is considered acceptable, given that the availability of transferable zeolite models that are appropriate as far as adsorption predictions are concerned for both polar and non-polar compounds is rather limited [16]. Furthermore, it is expected that the observed discrepancies between simulation and experimental results for this non-polar component are unable to cause significant errors in terms of mixture adsorption data, where the CH4 content of the gas phase is low and the adsorption of H2S is dominant. The calculated isosteric heat of adsorption data together with available experimental results for H2 and CH4 [5] are also plotted in Fig.1. The pressure- dependence of these data is weak. The general order of qH2S > qCH4 > qH2 is in line with expectations, bearing in mind that the isosteric heat of adsorption at low loadings proves the strength of interaction between the zeolite framework and the adsorbate molecules. q values are considerably higher for polar H2S than for non-polar substances and the order of magnitude of the former indicates the significance of electrostatic interactions. The relation of qCH4 > qH2 can be attributed to the greater polarizability of CH4 molecules (this is implicitly included in the attracting Lennard-Jones terms of the potential model), and to that the explicitly modelled real quadrupole moment of H2 molecules is very weak. Considering the greater uncertainties of these simulation results and that the experimental data for H2 and CH4 were obtained for zero coverage within a given temperature range, these simulation results also confirm the suitability of the models used in this study. Equilibrium adsorption selectivities were predicted for typical hydrodesulfurization stream outlets of petroleum refinery units separated by zeolite NaA at near-atmospheric pressures. The studied gas streams were comprised of between 1 and 2% H2S and ~95% H2; the remaining hydrocarbon content (low alkanes) was represented by the presence of CH4 in the model mixtures. For comparison, other compositions including very low and reasonably high H2S contents, as well as low pressure ranges were also investigated. The raw simulation results for the H2S-H2 mixtures in comparison with IAST predictions shown in Fig.2 illustrate well the dissimilar levels of adsorption of the two substances, with the exception of the nearly zero H2S contents of the bulk mixture. The IAST calculations underestimate the simulation results for H2 at higher pressures and on the whole overestimate the simulation results for H2S at lower pressures (for visual reasons, data obtained within the very low pressure range are not presented in this figure). The most accurate estimations were achieved at 10 kPa, which is an impractical parameter for the present applications. Strictly speaking, the hypothesis of IAST which states that the different adsorbate molecules have access to the same adsorbent surface cannot be applied to microporous adsorbents such as zeolite NaA, where the accessible surface area depends on the size of the adsorbate. In light of this, the IAST predictions can be considered to be remarkably accurate. a b c Figure 1. Equilibrium adsorption loading (n) and isosteric heat of adsorption (q) as a function of the bulk-gas pressure for pure hydrogen sulphide (a), hydrogen (b) and methane (c) on zeolite NaA at the temperatures indicated. The statistical uncertainties of the simulations results do not exceed the size of the symbols. The lines connecting simulation data at 298 K are only drawn to guide the eyes. sim.: simulation data; exp.: experimental data. KRISTÓF Hungarian Journal of Industry and Chemistry 12 The calculated equilibrium selectivities are defined as: b j b a j a nn nn S / / SH SH 2 2 , (4) where nj stands for the equilibrium number of moles of H2 in the investigated two-component mixtures or the sum of the equilibrium numbers of moles of H2 and CH4 in the three-component mixtures, as plotted in Fig.3. On the whole, this zeolite exhibits an exceptional level of selectivity of H2S to the other substances; this is not surprising given the significant differences between the equilibrium adsorption loadings of the pure components (cf. Fig.1). In the case of the two-component gas mixtures, the tendency of the data satisfies the criterion that at the low-pressure limit the selectivity as defined here should be independent of the composition of the bulk-gas mixture (it is the quotient of the ratio of the single-particle partition function of the two substances in the adsorbed phase and the ratio of their free-particle partition functions [34-35]). Because of technical reasons, at lower pressures the uncertainties of the selectivity data are relatively large as these data are calculated from simulation results at very low zeolite loadings. At higher pressures the separation efficiency of this zeolite is somewhat weaker. This and the fact that the change with pressure is less intense at lower H2S contents suggest the existence of a ‘crowding’ effect, which inhibits more strongly the sorption of the larger molecule, H2S. In the case of the investigated three-component mixtures, the overall picture is similar, but the selectivity values are smaller. This makes sense since the competitive effect of the additional component, CH4, for the adsorption sites is stronger. Yet, the values far in excess of 1000 obtained for the typical hydrodesulfurization streams (1-2% H2S and ~95% H2) are compelling. From the mixture adsorption data, once again it was verified that electrostatic effects control the adsorbent-adsorbate interactions with this zeolite, which implies that the amount of adsorption of pure H2 and CH4 should always be small. Adsorbed mole number data showed that the presence of the non-polar components does not affect the sorption of H2S in the adsorbed phase. This conclusion is also supported by the heat of adsorption data (not presented) calculated by assuming one-component mixtures (i.e. using Eq.(1)). These data turned out to be simply the amount-weighted average of the q values of pure components and are very close to qH2S at the given pressure. On the other hand, the degree of adsorption of the non-polar substances is reduced by the presence of H2S. Selectivities under real conditions (at 50, 100 and 200 kPa and with realistic H2S contents; 1, 2 and 5%) shown in Fig.4 make this fact obvious. Here, S values obtained from mixture adsorption simulations significantly exceed their counterparts calculated for an ideal case of independent adsorption (i.e. by substituting into Eq.(4) the pure- component adsorption loadings determined at pressures that are equal to the partial pressures of the mixture components). Extensive non-ideality in the adsorbed phase can also be seen from the comparative failure of IAST (which utilizes the assumption that the adsorbed mixture is an ideal solution) to predict the simulation results accurately at near-atmospheric pressures. It is remarkable that the selectivity of H2S to the two non-polar substances decreases as the temperature and partial pressure of H2S in the bulk gas increase. As the adsorption loading of the zeolite rises, steric hindrance plays an increasingly important role, and sorption of the larger H2S molecules reduces to a greater extent. The impact of an increase in temperature is as expected, e.g. from the higher qH2S values, but the magnitude of decrease in selectivity with temperature changes significantly as a function of the partial pressure of H2S. Data lines at the two investigated temperatures seem to converge to similar values at higher partial pressures, because the drop in the sorption of H2S as the temperature increases already becomes Figure 2. Comparison of the IAST predictions with simulation data for hydrogen sulphide (black) and hydrogen (blue). Equilibrium adsorption loading (n) as a function of the mole fraction of hydrogen sulphide (yH2S) in the binary gas phase mixture on zeolite NaA at 298 K and at the pressures indicated. The statistical uncertainties of the simulation results do not exceed the size of the symbols. (For interpretation of the references to colour in this figure, the reader is advised to refer to the online version of this article.) a b Figure 3. Equilibrium adsorption selectivity (S) on zeolite NaA at 298 K as a function of the bulk-gas pressure for binary (a) and ternary (b) gas mixtures with the compositions indicated. SELECTIVE REMOVAL OF HYDROGEN SULPHIDE USING ZEOLITE 45(1) pp. 9–15 (2017) 13 less significant at higher loadings. The two panels of Fig.4 also illustrate the above-mentioned difference between the data of the binary (a) and ternary (b) mixtures, i.e. the numerical values are smaller for the ternary mixtures. Besides this, the scatter of the points is larger in panel b, because the ratio of nH2 to nCH4 in the bulk gas unavoidably changes as the partial pressure of H2S increases (for compositions, see Fig.3). Finally, panel c in Fig.4 illustrates the influence of adsorbate- adsorbate attraction on the adsorption characteristics. It was simulated by eliminating the Coulomb potential and the attractive part of the Lennard-Jones potential, but retaining its soft-sphere repulsion potential component, when calculating the instantaneous adsorbate-adsorbate pair interactions in the adsorbed phase. The observed reduction in nH2S and selectivity is sizable enough to establish that like-like attraction is an important factor in the adsorption of H2S. 4. Conclusion In this work, molecular simulation predictions for the adsorption of H2S from simple non-polar gas mixtures of technological interest (hydrodesulfurization stream outlets in petroleum refinery units) on zeolite NaA were presented. The realistic all-atom intermolecular potential models adopted for the computations were validated by comparing the calculated isotherms of the pure substances with available experimental adsorption data. What is especially noticeable here is the matching of the experimental and simulated adsorption loadings for H2S that was achieved. The investigated zeolite exhibited a remarkable ability to capture H2S, from either binary or ternary mixtures with non-polar gases, namely H2 and CH4. The interactions between the polar H2S molecules and the hydrophilic zeolite framework were found to be particularly favourable, and the mixture-adsorption loadings for H2S essentially agreed with the corresponding pure component loadings (with the exception of the very low H2S contents of the bulk- gas mixture). The reverse is true when considering the adsorption of the non-polar gaseous components under technological conditions (at near-atmospheric pressures and with a small proportion of H2S in the bulk); their bindings to the inner surface sites of the zeolite were suppressed by H2S. These results can be of practical importance in terms of selectivity. Selectivities of H2S to the non-polar substances were generally higher at lower H2S partial pressures in the bulk gas, and well over 1000 for the range of H2S contents of the typical hydrodesulfurization streams. The obtained order of magnitude of the isosteric heat of adsorption data and the large decrease in selectivity with increasing temperature suggest that electrostatic interactions play a more pronounced role in the selective removal of H2S by zeolite NaA and the effect of size has only a limited impact. In association with this, it was also revealed that H2S-H2S attraction contributes to the preferred adsorption of this substance. Acknowledgement Present article was published in the frame of the project GINOP-2.3.2-15-2016-00053 (“Development of engine fuels with high hydrogen content in their molecular structures (contribution to sustainable mobility)”). We gratefully acknowledge the financial support of the Hungarian National Research Fund (OTKA K124353). The author would like to thank Tamás Kovács and Zoltán Ható (Department of Physical Chemistry, Institute of Chemistry, University of Pannonia) for their assistance in terms of data analysis. REFERENCES [1] Auerbach, S.M.; Carrado, K.A.; Dutta, P.K. (Eds.): Handbook of Zeolite Science and Technology (Marcel Dekker, New York) 2003 ISBN: 0-8247-4020- 3 a b c Figure 4. Equilibrium adsorption selectivity (S) as a function of the partial pressure of hydrogen sulphide in the bulk gas for selected binary (a) and ternary (b) gas mixtures on zeolite NaA at 298 K and 323 K, and at bulk gas pressures of 50 kPa, 100 kPa, and 200 kPa. Comparison of the mixture selectivity data with the selectivity data for independent adsorption (indep. ads.; panels a, b) and with mixture selectivity data calculated without adsorbate-adsorbate attractions (no attr., panel c). (For interpretation of the references to colour in this figure, the reader is referred to the web version of the article.) KRISTÓF Hungarian Journal of Industry and Chemistry 14 [2] Xu, X.; Yang, W.; Liu, J.; Chen, X.; Lin, L.; Stroh, N.; Brunner, H.: Synthesis and gas permeation properties of an NaA zeolite membrane, Chem. Commun. 2000 0, 603-604 DOI: 10.1039/B000478M [3] Aoki, K.; Kusakabe, K.; Morooka, S.: Separation of gases with an A-type zeolite membrane, Ind. Eng. Chem. 2000 39, 2245-2251 DOI: 10.1021/ie990902c [4] Okamoto, K.; Kita, H.; Horii, K.; Tanaka, K.; Kondo, M.: Zeolite NaA membrane: preparation, single-gas permeation, and pervaporation and va- por permeation of water/organic liquid mixtures, Ind. Eng. Chem. 2001 40, 163-175 DOI: 10.1021/ie0006007 [5] Zhu, W.; Gora, L.; van den Berg, A.W.C.; Kapteijn, F.; Jansen, J.C.; Moulijn, J.A.: Water va- pour separation from permanent gases by a zeolite- 4A membrane, J. Membrane Sci. 2005 253, 57-66 DOI: 10.1016/j.memsci.2004.12.039 [6] Yamamotoa, T.; Kimb, Y.H.; Kimb, B.C.; Endoa, A.; Thongprachana, N.; Ohmoria, T.: Adsorption characteristics of zeolites for dehydration of etha- nol: Evaluation of diffusivity of water in porous structure, Chem. Eng. J. 2012 181-2, 443-448 DOI: 10.1016/j.cej.2011.11.110 [7] Lee, S.H.; Moon, G.K.; Choi, S.G.; Kim, H.S.: Molecular dynamics simulation studies of zeolite- A. 3. Structure and dynamics of Na+ ions and water molecules in a rigid zeolite-A, J. Phys. Chem. 1994 98, 1561-1569 DOI: 10.1021/j100057a006 [8] Akten, E.D.; Siriwardane, R.; Sholl, D.S.: Monte Carlo Simulation of Single- and Binary Component Adsorption of CO2, N2, and H2 in Zeolite Na-4A, Energy & Fuels 2003 17, 977-983 DOI: 10.1021/ef0300038 [9] Furukawa, S.; Goda, K.; Zhang, Y.; Nitta, T.: Mo- lecular simulation study on adsorption and diffu- sion behavior of ethanol/water molecules in NaA zeolite crystal, J. Chem. Eng. Japan 2004 37, 67- 74 DOI: 10.1252/jcej.37.67 [10] Kristóf, T.; Csányi, É.; Rutkai, G.; Merényi, L.: Prediction of adsorption equilibria of water- methanol mixtures in zeolite NaA by molecular simulation, Mol. Sim. 2006 32, 869-875 DOI: 10.1080/08927020600934179 [11] Cosoli, P.; Ferrone, M.; Pricl, S.; Fermeglia, M.: Hydrogen sulphide removal from biogas by zeolite adsorption, Part I-II, Chem. Eng. J. 2008 145, 86- 92 and 93-99 DOI: 10.1016/j.cej.2008.07.034 & 10.1016/j.cej.2008.08.013 [12] Rutkai, G.; Csányi, É.; Kristóf, T.: Prediction of adsorption and separation of water-alcohol mix- tures with zeolite NaA, Microporous and Mesopo- rous Mat. 2008 114, 455-464 DOI: 10.1016/j.micromeso.2008.01.044 [13] Csányi, É.; Kristóf, T.; Lendvay, Gy.: Potential model development using quantum chemical in- formation for molecular simulation of adsorption equilibria of water-methanol (ethanol)mixtures in zeolite NaA-4, J. Phys. Chem. C 2009 113, 12225- 12235 DOI: 10.1021/jp902520p [14] Csányi, É.; Ható, Z.; Kristóf, T.: Molecular simula- tion of water removal from simple gases with zeo- lite NaA, J. Mol. Model. 2009 18, 2349-2356 DOI: 10.1007/s00894-011-1253-7 [15] Sun, Y.; Han, S.: Diffusion of N2, O2, H2S and SO2 in MFI and 4A zeolites by molecular dynamics simulations, Mol. Sim. 2015 41, 1095-1109 DOI: 10.1080/08927022.2014.945082 [16] Vujic, B.; Lyubartsev, A.P.: Transferable force- field for modelling of CO2, N2, O2 and Ar in all sil- ica and Na+ exchanged zeolites, Modelling Simul. Mater. Sci. Eng. 2016 24, 045002 (1-26) DOI: 10.1088/0965-0393/24/4/045002 [17] Pluth, J.J.; Smith, J.V.: Accurate redetermination of crystal structure of dehydrated zeolite A. Ab- sence of near zero coordination of sodium. Re- finement of Si, Al-ordered superstructure, J. Am. Chem. Soc. 1980 102, 4704-4708 DOI: 10.1021/ja00534a024 [18] Mikula, A.; Król, M.; Kolezynski, A.: Periodic model of an LTA framework, J. Mol. Model. 2015 21, 275 (1-9) DOI 10.1007/s00894-015-2820-0 [19] Bai, P.; Tsapatsis, M.; Siepmann, J.I.: TraPPE-zeo: Transferable Potentials for Phase Equilibria Force Field for All-Silica Zeolites, J. Phys. Chem. C 2013 117, 24375-24387 DOI: 10.1021/jp4074224 [20] Shah, M.S.; Tsapatsis, M.; Siepmann, J.I.: Devel- opment of the Transferable Potentials for Phase Equilibria Model for Hydrogen Sulfide, J. Phys. Chem. B 2015 119, 7041-7052 DOI: 10.1021/acs.jpcb.5b02536 [21] Kaminski, G.; Duffy, E.; Matsui, T.; Jorgensen, W.: Free Energies of Hydration and Pure Liquid Properties of Hydrocarbons from the OPLS All- Atom Model, J. Phys. Chem. 1994 98, 13077- 13081 DOI: 10.1021/j100100a043 [22] Darkrim, F.; Levesque, D.: Monte Carlo simula- tions of hydrogen adsorption in single-walled car- bon nanotubes, J. Chem. Phys. 1998 109, 4981- 4984 DOI: 10.1063/1.477109 [23] Waldman, M.; Hagler, A.: New combining rules for rare gas van der Waals parameters, J. Comput. Chem. 1993 14, 1077-1084 DOI: 10.1002/jcc.540140909 [24] Song, W.; Rossky, P.J.; Maroncelli, M.: Modelling alkane+perfluoroalkane interactions using all-atom potentials: Failure of the usual combining rules, J. Chem. Phys. 2003 119, 9145-9162 DOI: 10.1063/1.1610435 [25] Gubbins, K.E.; Quirke, N. (Eds.): Molecular Simu- lation and Industrial Applications: Methods, Ex- amples and Prospects (Gordon & Breach, Amster- dam) 1997 ISBN: 9056990055, 9789056990053 [26] Wolf, D.; Keblinski, P.; Phillpot, S.R.; Eggebrecht, J.: Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r-1 summation, J. Chem. Phys. 1999 110, 8254-8282 DOI: 10.1063/1.478738 [27] Demontis, P.; Spanu, S.; Suffritti, G.B.: Applica- tion of the Wolf method for the evaluation of Cou- lombic interactions to complex condensed matter systems: Aluminosilicates and water, J. Chem. Phys. 2001 114, 7980-7988 DOI: 10.1063/1.1364638 SELECTIVE REMOVAL OF HYDROGEN SULPHIDE USING ZEOLITE 45(1) pp. 9–15 (2017) 15 [28] Nicholson, D.; Parsonage, N.G.: Computer Simula- tion and the Statistical Mechanics of Adsorption (Academic Press, London) 1982 ISBN: 0125180608, 9780125180603 [29] Karavias, F.; Myers, A.L.: Isosteric heats of multi- component adsorption: thermodynamics and com- puter simulations, Langmuir 1991 7, 3118-3126 DOI: 10.1021/la00060a035 [30] Myers, A.; Prausnitz, J.M.: Thermodynamics of mixed-gas adsorption, AIChE J. 1965 11, 121-127 DOI: 10.1002/aic.690110125 [31] Simon, C.M.; Smit, B.; Haranczyk, M.: pyIAST: Ideal adsorbed solution theory (IAST) Python package, Comp. Phys. Commun. 2016 200, 364- 380 DOI: 10.1016/j.cpc.2015.11.016 [32] Cruz, A.J.; Pires, J.; Carvalho, A.P.; Carvalho, M.B.: Physical Adsorption of H2S Related to the Conservation of Works of Art: The Role of the Pore Structure at Low Relative Pressure, Adsorp- tion 2005 11, 569-576 DOI: 10.1007/s10450-005-5614-3 [33] Mohr, R.J.; Vorkapic, D.; Rao, M.B.; Sircar, S.: Pure and Binary Gas Adsorption Equilibria and Kinetics of Methane and Nitrogen on 4A Zeolite by Isotope Exchange Technique, Adsorption 1999 5, 145-158 DOI: 10.1023/A:1008917308002 [34] Tant, Z.; Gubbins, K.E.: Selective Adsorption of Simple Mixtures in Slit Pores: A Model of Me- thane-Ethane Mixtures in Carbon, J. Phys. Chem. 1992 96, 845-854 DOI: 10.1021/j100181a059 [35] Battisti, A.; Taioli, S.; Garberoglio, G.: Zeolitic imidazolate frameworks for separation of binary mixtures of CO2, CH4, N2 and H2: A computer simulation investigation, Microporous and Meso- porous Mat. 2011 143, 46-53 DOI: 10.1016/j.micromeso.2011.01.029 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 17–22 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0004 DYNAMIC SIMULATOR-BASED APC DESIGN FOR A NAPHTHA REDISTILLATION COLUMN LÁSZLÓ SZABÓ, 1 * KLÁRA KUBOVICSNÉ STOCZ, 1 LAURA SZABÓ, 1 SÁNDOR NÉMETH, 2 AND FERENC SZEIFERT 2 1 Department of Technology and Development, MOL Plc., Olajmunkás u 2., Százhalombatta, H-2443, HUNGARY 2 Department of Process Engineering, University of Pannonia, Egyetem u 10., Veszprém, H-8200, HUNGARY In this simulation study the operation of a naphtha redistillation column (a column with two feeds and three products) was analyzed with the application of Aspen HYSYS ® software. The simulator, structure of local con- trollers and the product-quality estimators were based on the data of an industrial column in the Danube Refin- ery. The aim of the analysis was to identify the dynamic and steady-state effects of heating and cooling as well as the sidestream of product qualities. The relationship between the tray temperatures and the quality of the products was also identified and inferential calculations were made. Based on the identified relationships, a two- level hierarchical control structure was developed. On the lowest level of the hierarchy are the local controllers of flowrates, liquid levels, pressure and duty. The inferential calculations are important components of the con- troller which serve as the controlled variables at the coordination level. The inputs of the estimators are the pro- cess data of the column, e.g. temperature, pressure and flowrate. On the top level of the control hierarchical structure the quality of the products are controlled by manipulating the setpoint of the local controllers. Based on the analysis of the Controlled Variable – Manipulated Variable relationship, closed-loop quality control was achieved with PID controllers. The result of the analysis may form the basis of Advanced Process Controller im- plementation. Keywords: distillation, multilevel control, MIMO, inferential, naphtha distillation, plant data 1. Introduction Distillation is the most widely used separation technique in the chemical and petrochemical industries. Around 95 % of the separation tasks of components are performed using distillation in the chemical industry, and the distillation units use approximately 3% of the total energy consumption of the world [1]. Therefore, the improvement of distillation controllers and processes is important because of their widespread use and huge levels of energy consumption. Distillation is a complex task because all of the equipment are Multiple Input, Multiple Output (MIMO) objects from a control point of view with strong interactions [2]. The performance of distillation controllers directly affects the productivity rates, utility usage and product quality. During the operation, attaining a sufficient level of product quality is the final goal. Hence, the precise calculation of the qualities is very important from an economic point of view. In this work, the multilevel control structure of the naphtha redistillation column has been analyzed (Fig.1). Furthermore, a simulator with the same local control *Correspondence: lszaboszhb@mol.hu structure as the existing distillation column in the Danube Refinery was created. The model was adapted according to the data of the sensors and laboratory measurements. The controlled quality values were provided by inferential calculations based on archived data of the unit. During the analyses a control structure was created which could be suitable for the quality control of the column, and also the resulting models can be used in the implementation of Advanced Process Controllers (APC). Figure 1. Structure of the control system SZABÓ, KUBOVICSNÉ STOCZ, SZABÓ, NÉMETH, AND SZEIFERT Hungarian Journal of Industry and Chemistry 18 2. Experimental For the analyses of dynamic responses a piece of commerical simulation software was used. Data evaluation and parameter identification of the inferential calculations were performed with R programming language. 2.1. Separation task The column has two feeds, the light naphtha feedstock possesses an initial boiling point of 12 °C and a final boiling point of 150 °C, while the boiling range of the heavy naphtha feedstock stretches from 57 to 197 °C. The task of the column is to separate the two feeds into three different naphtha products. The specification of the product streams are the following: distillate final boiling point should be between 58 and 65 °C, the final boiling point of the side product should be between 87 and 147 °C, and the initial boiling point of the bottom product should be between 99 and 103 °C. The simulator was implemented using Aspen HYSYS ® software. The structure of the simulator is shown in Fig.2. 2.2. Simulator The implementation of the column was achieved over two main steps. In the first step, a steady-state simulator of the column was developed. The structure of the simulator follows the process flow diagram in Fig.2. The column is divided into two parts, the light feedstock is loaded into the upper section (Column 107), and the heavy feedstock is fed into the lower section (Column 102). The towers, air coolers and heat exchangers were simulated with rigorous blocks and the model was adapted to data of the sensors and laboratory measurements. In the second step, the dynamic simulator was created based on the steady-state model. During the implementation of the dynamic simulator the calculation blocks were defined from the sizes of the actual equipment, e.g. the parameters of the trays in the tower. In some cases, pseudo pieces of equipment had to be defined to achieve a more realistic simulator, e.g. the bottom of the columns were modelled with a vertical drum to create a more detailed model. 2.3. Measurement Disturbances Before determining the parameters of the controller, noise and time delays were applied to the calculated values. The disturbaces were determined based on our previous experience. Transfer function blocks were used to simualte measuring instruments. Their parameters can be seen in Table 1. These signals were used later as controlled values. 2.4. Control Structure of Column Fig.1 represents the overall structure of controllers. On the lower level of the hierarchy the local controllers can be found. On this level, the controllers maintain the operating conditions (pressure and liquid level controllers) of the column and eliminate the disturbances (mass flowrate controllers) of the environment. The controller output variables are the valve positions and the process variables are the mass flowrates, liquid levels and pressure. The structure of the local controllers can be seen in Fig.2, which is the same as the structure of the controller in the real column. The FC1 and FC2 controllers eliminate the fluctuation of the feed mass flowrate. The PC and the TC controllers compensate for weather changes and the fluctuation of fuel gas. The liquid level controllers (LC1, LC2 and LC3) were necessary to ensure normal operation because the refluxes and the boil-up flowrate would be zero if there is no liquid phase in the reflux drum and at the bottom of the columns [3]. PI controllers were applied at the local level. The degree of separation is based on the sustenance of the temperature difference between Table 1. Disturbances Value type Noise Delay Lag Mass flowrate 1.0 % 15 s 0 s Pressure 0.2 % 0 s 0 s Temperature 0.3 % 30 s 60 s Level 1.0 % 60 s 0 s Figure 2. Structure of the simulator DYNAMIC SIMULATOR-BASED APC DESIGN 45(1) pp. 17–22 (2017) 19 products. The main target of the operation is to achieve specified product qualities. The provision of online quality values for the control is an important task in the industry. One of the feasible solutions is to apply inferential calculations. Generally the inputs of these calculations are temperatures, because the quality of products is strongly correlated with the tray temperatures of the columns. Consequently, the temperature and quality controls serve the same purpose. Therefore, in this paper the quality of the products was controlled at the top level of the control hierarchy. The controller outputs on this level were the setpoints of the local controllers (FC3, FC5 and TC) and the controlled values were the outputs of the inferential calculations. 3. Results and Analysis 3.1. Validation of the Simulator In the first analysis, the dynamic behavior of the simulator was validated. (The steady-state accuracy was validated during the construction of the model.) The model was compared to the historical data of the real column. As Fig.3 shows, the coil outlet temperature (COT) of the furnace was changed over two steps. During this process the top temperature of Column 107 was controlled. The head and bottom temperatures of Column 102 were analyzed, because the effect of the COT was significant on these parameters. The operating point of the real column was different from the simulated one, however, in this analysis only its dynamic behavior was evaluated. Therefore, in the simulator the setpoint of the COT was changed in the same direction and by the same amount as the measured data. The signals were normalized in order to compare the trends easily. It can be seen that the curves possess the same correlation in both the simulated and real columns. The dynamics of the responses were also identical in both cases, thus the simulator was found reliable for further analyses. 3.2. Analysis of Local Controllers For the tuning of the local controllers previous experiences were drawn on [4]. It was assumed that the controllers exhibit no significant interactions so they can be tuned independently. Before determining the parameters of the controllers, filters were applied, hence, the dynamics of the loops were approximated to the real control loops. The control loops were tuned in the following order: flowrate (FC1-7), the liquid level (LC1-3) and pressure (PC) controllers. In the last step the parameters of the coil outlet temperature controller (TC) were determined. For the loops, Honeywell equation B [5-6] was used, and the proportional and integrating parameters were determined using the built-in autotuner tool. After automatic tuning, the loops were tested and fine tuning was carried out to avoid fluctuation of the loops. The behavior of the FC1 is presented in Fig.4. As the figure shows, the controllers follow the setpoint without exceeding it, which ensures a robust local controller system. The raw calculated signal, the signal loaded with noise and the filtered (controlled) signals are also presented in Fig.4. As the results show, the control of the raw values is stable and is the first signal to reach a setpoint. 3.3. Inferential Calculations The quality estimators were created for the final boiling points of the distillate and the side product and for the initial boiling point of the bottom product. For the inferential fitting the laboratory and process data was used; the data was collected over a period of one year. In the case of the process data, all relevant measurement points of the column were considered (flowmeters, Figure 4. FC1 controller Figure 3. Validation of the Simulator SZABÓ, KUBOVICSNÉ STOCZ, SZABÓ, NÉMETH, AND SZEIFERT Hungarian Journal of Industry and Chemistry 20 temperature and pressure measurements). One hour- long periods were aggregated from the online data before each laboratory sample time. From these values pressure compensated temperatures (PCT) and ratios of the mass flowrates were also calculated. The inputs of the equations (see Table 2) were determined using correlation analysis, thus the qualities were calculated from the most relevant values. The parameters of the equations were calculated using linear regression. The results of the fittings are presented in Fig.5. 3.4. Structure of Quality Controllers According to the separation task the controlled variables (CV) are: distillate final boiling point (CV1), side- product final boiling point (CV2) and bottom product initial boiling point (CV3). The manipulated variables are the following: setpoints of FC3 (MV1), FC5 (MV2) and TC (MV3). The structure of the quality controllers was determined with the help of the relative gain array (RGA) [7]. To create the array, a step test was performed for each manipulated variable (MV) in the simulator. The RGA is presented in Table 3. The elements of the matrix evaluate a possible loop with the given MV – CV pair. If the element of the array is negative, closing other loops will reverse pairing gain. This configuration will become unstable if the other loops operate. If the element of the RGA exceeds one, interaction with other loops has an inhibiting effect. Closing other loops decreases pairing gain. If the value is between 0 and 1 then pairing gain will increase by closing the other loops; the interaction is the greatest at 0.5. In the case of this pairing the loop should be tuned when the other loops are automated. If the RGA value is 1 then the other loops have no effect on the pairing gain. In the case of the CV1 the evident manipulated variable is MV1. The other two can be paired in two ways. According to the RGA, the obvious pairing should be CV2 – MV2 / CV3 – MV3. However, there is a primary effect between MV2 and CV2, and it possesses a negative gain. This behavior could make the control loop instable, or the settling time of the controlled object would be large. The dynamic responses are presented in Fig.6. In this analysis the MVs were stepped during the first minute. Only the local controllers were in closed loops during the analyzed period. Table 2. Input parameters of inferential calculations Inferential Input 1 Input 2 Final boiling point of distillate PCT of 107 head temperature Final boiling point of side product PCT of 107 bottom temperature ���ℎ� ���ℎ�ℎ� ���ℎ� ���� Initial boiling point of bottom product PCT of 102 bottom temperature Figure 6. Open loop test CV2 – MV2 Table 3. Relative Gain Array MV1 MV2 MV3 CV1 7.7 -1.2 -5.5 CV2 -6.9 1.9 6.0 CV3 0.2 0.3 0.5 Figure 5. Performance of inferential calculations DYNAMIC SIMULATOR-BASED APC DESIGN 45(1) pp. 17–22 (2017) 21 The following final control pairings were applied: CV1 – MV1, CV2 – MV3 and CV3 – MV2. According to this analysis quality control can be achieved using simple PID control loops. Therefore, the analysis of the closed-loop quality control in the following chapter was performed using only PID loops. This analysis was also the basis of an APC implementation. The initial APC models can be seen in Fig.7. 3.5. Operation of Quality Controllers In this analysis the operation of the quality controllers was analyzed in the case of a disturbance. The mass flowrate of the light naphtha feed was changed at 5 o’clock. During the analyzed period all of the local and quality controllers were in closed loops. According to the results in Fig.8, the controller system can compensate for this change. In the figure the values of the MVs were normalized between their minimum and maximum limits. In the case of CV1 this quality did not violate the specification limit, and the controlled variable was close to the setpoint after 2.5 hours. This parameter possessed the smallest peak of the three CVs at approximately 1.3 °C. CV2 exhibited a peak at approximately 2 °C, nevertheless, this was well within the specification limit so it did not violate it. The settling time was about 7.5 hours. CV3 exhibited a peak at approximately 3 °C which violated the specification limit. After 2.5 hours the CV once again exceeded the limit. The settling time was ~15 hours. As can be seen the outputs (OP) were moved after the settling of the CVs, but because of the long settling time the whole time period was not presented. Figure 8. Operation of Quality Controllers Figure 7. APC models SZABÓ, KUBOVICSNÉ STOCZ, SZABÓ, NÉMETH, AND SZEIFERT Hungarian Journal of Industry and Chemistry 22 4. Conclusion For the analysis of the naphtha redistillation column at MOL Plc. the dynamic simulator of the unit was adapted to the data of the real equipment. The steady- state accuracy and the dynamic behavior of the model were evaluated. The simulator is reliable for the analysis of the column. For the operation of the distillation column a two- level control structure was implemented in the simulator. On the lower level of the control hierarchy controllers were created which facilitated normal operating conditions of the column (pressure and liquid- level controllers) and eliminated the disturbances of the environment (mass flowrate controllers). To ensure quality signals, inferential calculations were created based on the historical data of the real column. The inputs of the calculations were mainly the temperature measurements of the columns. Hence, the controlled temperatures were chosen during the correlation analysis of the inferential calculations; on the top level of the control hierarchy the controlled variables were the outputs of the quality estimators. The manipulated variables were the setpoints of the local controllers. The pairing of the manipulated and controlled variables should be defined carefully. The first aspect of the pairing was the gain between the variables, however, the dynamic responses were also considered. To eliminate the interaction and to improve the performance of the quality controller, decoupling [8] or a model predictive controller [9] was recommended. The simulator-based column analysis can be a good foundation for an Advanced Process Controller. The CV-MV pairing and resulting single input-single output models can be used as initial models in Advanced Process Controllers. Furthermore, the simulator can be used for the analysis of column operation issues in the plant. The resulting inferential calculations have already been implemented and are available for the operators to facilitate quality control. Acknowledgements We acknowledge the financial support of Széchenyi 2020 under the project EFOP-3.6.1-16-2016-00015. REFERENCES [1] Mansour, K.; Mansour, E.: Rigorous optimization of heat-integrated and Petlyuk column distillation configurations based on feed conditions, Clean Techn. Environ. Policy, 2009 11(1), 107–113 DOI: 10.1007/s10098-008-0171-6 [2] Skogestad, S.; Morari, M.: The dominate time con- stant for distillation columns, Computers and Chemical Engineering, 1987 11(6), 607–617 DOI: 10.1016/0098-1354(87)87006-0 [3] Buckley, P.S.: Techniques of Process Control, (Wiley, London, United Kingdom) 1964 ISBN- 10:0882757776 [4] Szabó, L.; Németh, S.; Szeifert, F.: Three-Level Control of a Distillation Column, Engineering, 2012 4(10), 675–681 DOI: 10.4236/eng.2012.410086 [5] Honeywell: Application Module Algorithm Engi- neering Data, 1-800 343-0228 [6] Aspentech: Aspen HYSYS Unit Operations Refer- ence Guide, Version Number: V10 [7] Skogestad, S.: Dynamics and Control of Distilla- tion Columns: A Tutorial Introduction, Chemical Engineering Research and Design, 1997 75(6), 539–562 DOI: 10.1205/026387697524092 [8] Benz, S.J.; Scennal, N.J.: An Extensive Analysis on the Start-up of a Simple Distillation Column with Multiple Steady States, The Canadian Journal of Chemical Engineering, 2002 80(5), 865–881 DOI: 10.1002/cjce.5450800510 [9] Luyben, W.L.: Process Modeling, Simulation, and Control for Chemical Engineers (McGraw-Hill, New York, USA) 1990 ISBN-10:0070391599 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 23–27 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0005 REPLACEMENT OF BIASED ESTIMATORS WITH UNBIASED ONES IN THE CASE OF STUDENT'S t-DISTRIBUTION AND GEARY’S KURTOSIS GERGELY TÓTH* AND PÁL SZEPESVÁRY Institute of Chemistry, Faculty of Science, Eötvös Loránd University, 1/A Pázmány Péter sétány, Budapest, H-1117, HUNGARY The use of biased estimators can be found in some historically and up to now important tools in statistical data analysis. In this paper their replacement with unbiased estimators at least in the case of the estimator of the population standard deviation for normal distributions is proposed. By removing the incoherence from the Stu- dent’s t-distribution caused by the biased estimator, a corrected t-distribution may be defined. Although the quantitative results in most data analysis applications are identical for both the original and corrected t- distributions, the use of this last t-distribution is suggested because of its theoretical consistency. Moreover, the frequent qualitative discussion of the t-distribution has come under much criticism, because it concerns artefacts of the biased estimator. In the case of Geary’s kurtosis the same correction results (2/π)1/2 unbiased estimation of kurtosis for normally distributed data that is independent of the size of the sample. It is believed that by re- moving the sample-size-dependent biased feature, the applicability domain can be expanded to include small sample sizes for some normality tests. Keywords: unbiased estimator, normal distribution, Anscombe-Glynn test, Jarque-Bera test, Bonett-Seier test 1. Introduction The Student's t-distribution [1] is one of the most widely used statistical functions in experimental practice. It is well documented that experts active in analytical, physi- cal and clinical chemistry, biology, agriculture, ecology, economy as well as forensic science or even legal repre- sentatives apply this tool to formulate solid statements, e.g. population means, to differentiate between two sample means, etc., which are in most cases necessarily based on a limited number of observations. Fortunately, scientists are in this respect soundly supported by a lot of textbooks, standards, software, software systems, and last but not least trained in elements of statistics. How- ever, it should be a moral obligation to be aware of the evolution of this routinely used method, its principles and often tacitly supposed assumptions. Although not crucial in daily use, it is worthwhile to know that be- sides Student’s t-distribution there are different func- tions which may be suitable for the determination of percentiles in the same way as the t-distribution. They may differ mainly in terms of alternate estimators for the population mean and population standard deviation. An attractive variant is presented in this work emphasiz- ing its theoretically fully consistent feature on the con- trary to the Student's t-distribution. *Correspondence: toth@chem.elte.hu The Student’s t-distribution corresponds to a ratio of normally distributed random variables to chi- distributed random variables (see the references in the historical review of Zabell [2]). Chi-distributed varia- bles postulate normally distributed data as well. Gosset [1] used in his definition an estimate of the standard deviation which is biased relating the population stand- ard deviation and even the variance in the case of nor- mally distributed random variables as was shown earlier by Helmert [3-4]. Amazingly, when Fischer proposed a transformation of Gosset’s original z variable to 1 Nzt [5-6], he chose an estimate for the standard deviation which is also biased relating σ. A corrected tc- distribution is proposed that fulfils all theoretical re- quirements and yields a more normal distribution-like shape. It is consequently based on normal sample data and uses an unbiased estimator. A similar correction can be applied to Geary’s kur- tosis that is the ratio of the mean deviation to the stand- ard deviation [7]. In this work the use of the correction is proposed in order to eliminate the sample-size de- pendency of the mean Geary’s kurtosis on normally distributed data. Finally, some remarks are made on statistical tests based on sample-size dependent values in order to ex- tend their applicability to small sample sizes. TÓTH AND SZEPESVÁRY Hungarian Journal of Industry and Chemistry 24 2. Theoretical Background The square root of the mean of the square of the centred observations of the random variable, Y,   N yy s N i i N     1 2 (1) may be an estimator of the population standard devia- tion, σ, of a sample. N is the sample size and y denotes the sample mean. Bessel pointed long before to the bias of sN and proposed:   1 1 2      N yy s N i i . (2) Although �� is already an unbiased estimator of the variance, ��, irrespective of the distribution of Y, the statistics �� and � are both biased in terms of the stand- ard deviation, �. However, there is a correction [5-6, 8- 9]: , 2 1 2 1 2 )(4                  N N N Nc (3) which applies to � by transforming it, in the case of normally distributed variables, into an unbiased statistic:  Nc s s 4 c  . (4)  Nc4 follows from Helmert’s papers or Cochran's theorem [10]. For normally distributed Y, )(1 sN  exhibits a chi distribution with N-1 degrees of freedom.  Nc4 is the expected value of s/σ. There is a  Nc2 value as well, that is equal to the expected value of �� �⁄ . The correcting effect of  Nc4 may be consider- able for low values of N (Table 1). Table 1. Selected c4(N) corrections N c4(N) 2 0.7979 3 0.8862 4 0.9213 5 0.9400 10 0.9727 20 0.9869 30 0.9914 Figure 1. Student’s t-distribution compared with the standard normal distribution The correction term is frequently used in Statistical Process Control (SPC) to define ±3σ intervals and pro- cess standard deviations that are determined for samples of different sizes. 3. Results and Discussion 3.1. Discussion of the density functions of the Student’s t-distribution Let Y be an independent N(μ,σ2) variable. The Ns y t /   (5) statistic exhibits the usual Student t-distribution with N- 1 degrees of freedom. The density function (Eq.(6)) can be derived as the quotient of normal and chi-distributed random variables. The independency of the nominator and denominator can be shown and it is also proved by a series of theorems that t defined in Eq.(5) follows a Student’s t-distribution with N-1 degrees of freedom (see the references in [2]):   22 1 1 2 1 2 1 1 )( N N t N N N tf                              (6) The distribution sketched in Fig.1 is widely known and needs no comments except for its flaw: it is based on the biased statistic (Eq.(2)) for σ. By replacing the standard deviation of the sample, s, in Eq.(5) with the unbiased equivalent one, sc, cor- rected by c4(N), results in a new value: Ns y t /c c   , (7) and in the density function: NOTE ON STUDENT’S T-DISTRIBUTION 45(1) pp. 23–27 (2017) 25       2 2 4 2 c 2 2 4 2 c 4 c 1 1 2 1 1 1 2 1 1 2 )( 1 )( N N Nc t Nc t N N N Nc tf                                         . )( )( )1,2/(crit 4 )1,2/(crit 4 )1,2/(crit cc     N NN st Nct Nc s ts   Table 2. Sample-size dependence of Geary’s and Pearson’s kurtosis. The c4(N)-corrected Geary’s and the size-corrected Pearson’s kurtosis [11] provided the ∞ limit values for all sample sizes within the statistical uncertainty of the 105 random samples. N Geary’s kurtosis Pearson’s kur- tosis 3 0.9004 1.5000 4 0.8659 1.8005 5 0.8489 1.9996 6 0.8385 2.1437 7 0.8319 2.2501 8 0.8267 2.3340 9 0.8233 2.4006 10 0.8200 2.4557 20 0.8085 2.7141 50 0.8021 2.8820 100 0.7999 2.9408 ∞ 0.7979 3.0000 (8) which is shown in Fig.2. As expected, the corrected Student’s t-distribution (Student’s tc-distribution) con- sists of more distinct peaks, it exhibits fatter tails than the standard normal distribution, but rather oddly at its maximum, i.e. when tc = 0, the value f(tc) does not de- pend on N and is represented by the value 2/1 , the same as for those of normal distribution. The function f(tc) and the Student’s t-distribution without a doubt differ in this respect. The obvious differences between the Student and modified Student distributions do not complicate their daily usage. The confidence intervals calculated using s and t or sc and tc estimators and distributions, respectively, N ts y N )1(,2/ cc    and N st y N )1(,2/    (9) do not differ, because (10) Evidently, corresponding estimators and functions should be used. 3.2. The effect of the correction to Geary’s kurtosis Geary’s kurtosis [7], wN, is the ratio of the mean abso- lute deviation (MAD) to the standard deviation (Eq.(11)). It is an alternative to Pearson’s kurtosis based on the fourth moment. The expected value of Geary’s kurtosis depends on the sample size even for normally distributed data [7]. The mean values of 105 random samples from standard normal distributions are shown in Table 2. N N i i N s yy N w     1 1 (11) If the c4(N) corrections (Table 1) are applied during the calculation of the nominator of the ratio as wN,corr=wN c4(N), the expected value of the kurtosis is (2/π)1/2 for all sample sizes. This means that the platy- kurtic and the leptokurtic features of a sample can be found without searching for the size-dependent dividing value in tables. 3.3. Sample-size bias in statistical tests Geary’s kurtosis and its transformed values are used in normality tests due to their enhanced sensitivity to some leptokurtic deviations from normality [12]. Contrary to the case of the Student’s t-distribution, where the cor- rection has no effect on the t-test, here the effect of the correction is not cancelled. Generally, the size depend- ence decreases the performance of the tests for small sample sizes. This feature is interpreted by users as a recommendation that the tests are unsuitable for small sample sizes. In the same way, neglect of sample-size dependence is applicable in tests where Pearson’s kurto- sis is used. The calculated mean value of Pearson’s kur- tosis is shown in Table 2 but its convergence is rather weak to the theoretical value of 3. It should be noted here that the sample-size unbiased estimator of kurtosis can be easily calculated [11]. Figure 2. Corrected Student’s t-distribution compared with the standard normal distribution TÓTH AND SZEPESVÁRY Hungarian Journal of Industry and Chemistry 26 Table 3. Sample-size dependence of five normality tests based on unbiased or sample-size-dependent biased estimators. The numbers show the ratio of the rejected null hypotheses to all trials at a significance level of 0.05 from 105 random samples. N Shapiro- Wilk D’ Agostino Anscombe- Glynn Bonett- Seier Jarque -Bera 4 0.0502 5 0.0521 0.0373 6 0.0477 0.0189 0.0087 7 0.0492 0.0343 0.0282 8 0.0505 0.0365 0.0348 9 0.0507 0.0538 0.0369 0.0390 0.0024 10 0.0498 0.0528 0.0394 0.0417 0.0058 20 0.0497 0.0525 0.0406 0.0421 0.0089 50 0.0498 0.0500 0.0466 0.0470 0.0241 100 0.0499 0.0493 0.0533 0.0490 0.0368 In Table 3 the type-I error of some normality tests calculated on 105 standard normal samples is shown. In the calculation the ‘moments’ package in R was used [13]. Table 3 contains the ratio of the samples to all samples where the H0 hypothesis of normality was re- jected at the significance level of 0.05. The Shapiro- Wilk method [14] uses the ratio of two unbiased estima- tors of the variance, and is suitable for all data sizes. The skewness test of D’Agostino [15] slightly overesti- mates the number of rejected cases. It is based on the normalized third-order central moment definition of skewness, where the expected value is estimated with- out bias. The Anscombe-Glynn [16] test applies Pear- son’s kurtosis without size correction using a biased estimation of normally distributed data. The Bonett- Seier [12] test shown here uses Geary’s kurtosis and the Jarque-Bera test [17] combines skewness and Pearson’s kurtosis. The performance of the last three tests was rather weak for small sample sizes, whereof one cause might be the lack of correction for small sample sizes even for normally distributed data. These tests are usu- ally only recommended for medium and large sample sizes. The correction should extend the applicability domain to small sample sizes. Of course, the type-I er- ror of normally distributed data is only one narrow as- pect of a test, detailed analysis should be performed to investigate the effect of correction on many distribu- tions, like, e.g. in ref. [12]. 4. Conclusion Nowadays, data are evaluated by computers and biased estimators can be replaced by unbiased ones, even if their calculation schemes are complex. It has been shown that, in terms of Student’s t- distributions, to decide upon the confidences of statis- tics one has two functions which are completely equiva- lent as far as practical applicability is concerned. They can, however, be distinguished theoretically. The asser- tion that only the unbiased estimator should be recog- nized as the correct one implies the use of the corrected Student’s t-distribution, f(tc). In that case the known shape of the Student’s t-distribution may be labelled as an artefact and the usual application of the Student’s t- distribution as a production of “correct numbers by an incoherent theory”. In the case of Geary’s kurtosis, the correction re- moves the sample-size dependence from the expected value. This change of distinguishing platykurtic or lep- tokurtic features of the sample is simpler than using the original version of Geary’s kurtosis. Furthermore, sub- tracting (2/π)1/2 results in a number to be interpreted in a similar way to the excess kurtosis obtained by subtract- ing 3 from the Pearson’s kurtosis. As a further study, the use of unbiased/sample- size-dependent corrections to extend the applicability domain to small sample sizes in the case of normality tests is recommended. It is believed that the use of bi- ased estimators was acceptable before the age of com- puters and a systematic change to unbiased ones might be necessary in terms of statistics and standards with regard to industrial processes. Acknowledgement The authors are sincerely grateful for the valuable comments of Prof. S. Kemény and other participants during the scientific discussion after the presentation concerning Student’s t-distribution at the SCAC 2010 conference in Budapest in September 2010. REFERENCES [1] Student (Gosset, W.S.): The Probable Error of a Mean, Biometrika 1908 6(1), 1–25 DOI: 10.2307/2331554 [2] Zabell, S.L.: On Student’s 1908 Article “The Prob- able Error of a Mean”, J. Am. Stat. Assoc. 2008 103(481), 1–7 DOI: 10.1198/016214508000000049 [3] Helmert, F.R.: Über die Wahrscheinlichkeit der Potenzsummen der Beobachtungsfehler und über einige damit in Zusammenhang stehende Fragen, Z. Math. Phys. 1876 21, 192–218 [4] Helmert, F.R.: Die Genauigkeit der Formel von Peters zur Berechnung des wahrscheinlichen Beo- bachtungsfehles directer Beobachtungen gleicher Genauigkeit, Astronomische Nachrichten 1876 88(8-9), 113–131 DOI: 10.1002/asna.18760880802 [5] Fischer, R.A.: Applications of “Student’s” distribu- tion, Metron 1925 5, 90–104 [6] Fischer, R.A.: Statistical methods for research workers, (Oliver & Boyd, Edinburgh and London) 1925. [7] Geary, R.C.: Moments of the ratio of the mean de- viation to the standard deviation for normal sam- ples, Biometrika 1936 28(3/4), 295–307 DOI: 10.2307/2333953 [8] Grubbs, F.E.; Weaver, C.L.: The best unbiased es- timate of population standard deviation based on group ranges, J. Am. Stat. Assoc. 1947 42(238), 224–241 DOI: 10.1080/01621459.1947.10501922 NOTE ON STUDENT’S T-DISTRIBUTION 45(1) pp. 23–27 (2017) 27 [9] Vincze, I.: Matematikai statisztika ipari alkal- mazásokkal, (Műszaki Kiadó, Budapest), 1975 pp. 58, 89, 165–168 ISBN 963-10-0472-4 [10] Cochran, W.G.: The distribution of quadratic forms in a normal system, with applications to the analy- sis of covariance, Math. Proc. Cambridge 1934 30(2), 178–191 DOI: 10.1017/S0305004100016595 [11] Joanes, D.N.; Gill, C.A.: Comparing measures of sample skewness and kurtosis, J. Roy. Stat. Soc. D- Sta. 1998 47(1), 183–189 DOI:10.1111/1467-9884.00122 [12] Bonett, D.G.; Seier, E.: A test of normality with high uniform power, Comput. Stat. Data An. 2002 40(3), 435–445 DOI: 10.1016/S0167-9473(02)00074-9 [13] Komsta, L.; Novomestky, F.: Moments R package, version 0.14, 2015 at http://cran.r-project.org, ac- cessed in October 2017 [14] Shapiro, S.S.; Wilk, M.B.: An Analysis of Vari- ance Test for Normality (Complete Samples), Bio- metrika 1965 52(3/4), 591–611 DOI: 10.1093/biomet/52.3-4.591 [15] D’Agostino, R.B.: Transformation to Normality of the Null Distribution of G1, Biometrika 1970 57(3), 679–681 DOI: 10.1093/biomet/57.3.679 [16] Anscombe, F.J.; Glynn, W.J.: Distribution of the kurtosis statistic b2 for normal samples, Biometrika 1983 70(1), 227–234 DOI: 10.1093/biomet/70.1.227 [17] Jarque, C.M.; Bera, A.K.: Efficient tests for nor- mality, homoscedasticity and serial independence of regression residuals, Economic Letters 1980 6(3), 255–259 DOI: 10.1016/0165-1765(80)90024-5 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 29–36 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0006 FORMATION, PHOTOPHYSICS, PHOTOCHEMISTRY AND QUANTUM CHEMISTRY OF THE OUT-OF-PLANE METALLOPORPHYRINS ZSOLT VALICSEK, 1* MELITTA P. KISS, 1 MELINDA A. FODOR, 1 MUHAMMAD IMRAN, 2 AND OTTÓ HORVÁTH 1 1 Department of General and Inorganic Chemistry, Institute of Chemistry, Faculty of Engineering, University of Pannonia, Egyetem u. 10., H-8200 Veszprém, HUNGARY 2 Department of Chemistry, Baghdad-ul-Jadeed Campus, The Islamia University of Bahawalpur, 63100 Bahawalpur, PAKISTAN Among the complexes of porphyrins, special attention has been paid to those possessing out-of-plane (OOP) structures, for the formation of which the size, as well as the coordinative character of the metal center are re- sponsible. In these coordination compounds, the central atom cannot fit coplanarly into the cavity of the ligand, hence, it is located above the porphyrin plane, distorting it. Equilibria and kinetics of the complex formation, spectrophotometric, photophysical and primary photochemical properties of post-transition and lanthanide OOP metalloporphyrins were investigated, in addition electronic structural calculations were performed; hence, the general OOP characteristics were determined.Meanwhile, few doubtful questions have attempted to be an- swered concerning the categorization of metalloporphyrins, the borderline case complexes and hyperporphyrins. Keywords: out-of-plane metalloporphyrins, formation kinetics, UV-Vis spectrophotometry, photo- chemistry, borderline case complexes 1. Introduction Porphyrins and their derivatives play important roles in several biochemical systems. Four pyrroles are connected through methylidine bridges, forming the porphin ring. Its planar structure with an extended conjugated π-electron system provides aromatic characteristics and a special coordination cavity for the binding of metal ions of suitable radius [1-2]. Metalloporphyrins are the central parts of naturally important compounds, e.g., magnesium(II) chlorins in bacteriochlorophylls and chlorophylls; iron(II) protoporphyrin in hemoglobin; and iron(III) protoporphyrins in myoglobin, cytochromes, oxidase, peroxidase, catalase, and oxoanion reductase enzymes. Ringed tetrapyrroles provide strong chelating effect which can promote the hyperaccumulation of rare metal ions in living cells, and also in abiotic environments, e.g. in kerogens. In porphyrins the conjugation favors a planar structure. However, peripheral substituents or the metal center (originating from its size or axial ligand) can cause geometrical distortion. This certainly affects the functions of enzymes, as well as the biosynthesis of metalloporphyrins. In chemical research, due to distortion, redox potentials, basicity, reactivity, catalytic *Correspondence: valicsek@almos.uni-pannon.hu activity and coordinative abilities of porphyrins can be modified. Also, due to the distortion, the degree of symmetry decreases, resulting in characteristic spectral changes in various ranges of the electromagnetic spectrum. The most frequent types of distortions are dome, saddle, ruffled and wave (chair-like) [3]. Overcrowded substitution on the periphery [3-4] or insufficiently short metal-nitrogen bonds due to the shrinkage of the coordination cavity can cause the ruffled or saddled deformation [2, 5-7]. If, however, the M-N bonds are significantly longer than half the length of the diagonal N-N distance in the free-base porphyrin, dome deformation can occur. This happens if the radius of the metal center exceeds the critical value of about 75-90 pm (depending on the type of porphyrin ligand) or square planar coordination is not preferred. Such metal ions are too big to fit into the ligand cavity. Hence they are located above the plane of the pyrrolic nitrogens; forming sitting-atop (SAT) or out-of-plane (OOP - see Fig.1) complexes, displaying thermodynamic instability, kinetic lability, typical photophysical features and photochemical reactivity [8- 9]. In this work, we review our recent results regarding the formation, structure and photoinduced behavior of water-soluble OOP metalloporphyrins. These complexes with a diverse range of metal ions can be more simply produced in aqueous systems than in organic solvents. In this regard one of the most suitable free-base ligands is the anionic 5,10,15,20-tetrakis(4- sulfonatophenyl)porphyrin (H2TSPP4– - see Fig.1) due VALICSEK, KISS, FODOR, IMRAN, AND HORVÁTH Hungarian Journal of Industry and Chemistry 30 to its negative charge promoting the coordination of positively charged metal ions. Besides, this ligand is the most frequently used reagent among the free-base porphyrins [1]. 2. Experimental Analytical grade tetrasodium 5,10,15,20-tetrakis(4- sulfonatophenyl)porphyrin (C44H26N4O12S4Na4·12H2O = Na4H2TSPP·12H2O) (Sigma–Aldrich) and simple metal salts such as nitrate, sulfate, chloride or perchlorate were used for the experiments. The solvent was double-distilled water purified with a Millipore Milli-Q system. The pH of the majority of the metalloporphyrin solutions was adjusted to 8 by application of a borate buffer, whilst maintaining the ionic strength at a constant value of 0.01 M. In a few cases, the pH was regulated to 6, and the ionic strength to 1 M, by an acetate buffer, to hinder hydrolysis. The absorption spectra were recorded and the spectrophotometric titrations were monitored by using a Specord S-100 or a Specord S-600 diode array spectrophotometer. For the measurement of fluorescence spectra, a Perkin Elmer LS-50B or a Horiba Jobin Yvon FluoroMax-4 spectrofluorometer was applied. The latter piece of equipment supplemented with a time-correlated single photon counting (TCSPC) accessory was utilized to determine fluorescence lifetimes, too. UV-Vis spectrophotometric data (molar absorption, fluorescence quantum yields and lifetimes) of the free-base porphyrin were used as references for the determination of those of metalloporphyrin complexes [1]. For the determination of photochemical properties via continuous irradiations, a piece of AMKO LTI photolysis equipment (containing a 200W Xe–Hg lamp and a monochromator) was applied. For the electronic structural calculations, the B3LYP Density Functional Theory (DFT) method and the LANL2DZ basis set were used. On the basis of our earlier quantum chemistry experiences, the sulfonatophenyl substituents exhibit negligible effects on the coordination of the metal center in the cavity; thus, the anionic porphyrin (H2TSPP4−) can be modeled on the unsubstituted porphin (H2P) [4, 10]. 3. Results and discussion 3.1. UV-Vis spectrophotometry Porphyrins and their derivatives belong to the strongest light-absorbing materials (both natural and artificial), therefore, ultraviolet-visible spectrophotometry is one of the most fundamental, in addition, most informative spectroscopic techniques in porphyrin chemistry. They possess two ππ* electronic transitions in the visible range of the electromagnetic spectrum: B- or Soret band at about 350-500 nm, usually with a molar absorbance of 105 M-1cm-1 (Fig.2), and Q bands at 500-750 nm generally with intensities of one order of magnitude less. These latter bands in free-base ligands split due to the presence of protons on two diagonally situated pyrrolic nitrogens, to be more precise, as a result of the reduced symmetry (because of the disappearance of the four-fold rotation axis) compared to the metallated or deprotonated forms. This split is not detectable in the Soret range, hence, these two types of bands in the visible region are remarkably different [1]. In the Soret region, compared to the corresponding free-base ligands, the typical in-plane metalloporphyrins (e.g. Fe3+, Au3+, Cu2+, Pd2+) exhibit blueshifts because the atomic orbitals of their metal centers which are covalently bonded in the plane can overlap more strongly with the highest occupied molecular orbitals (HOMO) of the ligand, resulting in a stronger reduction in energy; whereas the lowest unoccupied molecular orbitals (LUMO) do not change. Thus, the energy gaps between the excited and ground states become greater. In the OOP complexes, the atomic orbitals of the more weakly bonded metal ions (e.g. Cd2+, Hg2+, Tl3+) may slightly affect the unoccupied MOs and to a lesser extent the occupied ones, leading to a reduction of the energy gaps, i.e. an increase in the corresponding wavelengths (Scheme 1) [1, 9]. a) b) Figure 1. Structure of an in-plane metalloporphyrin {MTSPP=metallo-5,10,15,20-tetrakis(4- sulfonatophenyl)porphyrin} (a); and that of an out-of-plane complex (b) [9]. STUDY OF THE OUT-OF-PLANE METALLOPORPHYRINS 45(1) pp. 29–36 (2017) 31 Beside electronic factors, due to the rigidity of the porphyrins’ ringed structure, steric effects also influence the spectra: the redshift of absorption bands is one of the most common spectroscopic consequences of the non-planarity of porphyrin [3]. Octabrominated free- base porphyrin, H2TSPPBr8 4–, was applied to investigate the spectrophotometric effects of the macrocycle’s highly saddle-distorted structure (Fig.2). In porphyrins with aryl substituents, this distortion can lead to the extension of delocalization by the twisting of aryl substituents from a nearly perpendicular orientation closer to the porphyrin plane (Fig.1) [4]. The larger, post-transition metal ions, e.g. thallium(I), lead(II) and bismuth(III) ions, can cause a similarly large redshift of the porphyrins’ absorption bands. Since their complexes possess the most highly dome-distorted structures, also a ruffled-like deformation of the periphery superposes on this high degree of doming. Considering the spectral effects (bathochromic or not quite exactly hyperchromic effects), the complexes possessing such highly redshifted absorption bands used to be referred to as hyperporphyrins; depending on the highest occupied electron subshell of the metal center, p- or d-type hyperporphyrins. Previously in terms of this categorization of metalloporphyrins, only the electronic effects of the metal ion (through its electron configuration) were taken into consideration and not steric (distorting) effects [1]. Nevertheless, in the typical d-type hyperporphyrins, e.g. the low-spin chromium(III), manganese(III), nickel(II) and cobalt(III) porphyrins, the radius of the metal center, and thus, the metal-nitrogen bonds are too short, resulting in the contraction of the coordination cavity, along with the ruffled deformation of the macrocycle [2, 5-7]. 3.2. Equilibrium and kinetics of complex formation Porphyrins are peculiar ligands in terms of complexation due to their planar, cyclic, rigid, aromatic, tetradentate, as well as protonated structure. The formation of an out-of-plane complex of a large metal ion is usually at least two orders of magnitude faster than that of an in-plane one since a smaller metal ion is not able to coordinate to all four pyrrolic nitrogens of the reaction intermediate, in the cavity of which the two protons also remain {H2-P-M}. Therefore, dissociation of the metal ion is more favorable than that of the protons. Besides the insertion of a smaller metal ion, its dissociation from the in-plane complex of the end- product may be kinetically hindered due to the rigidity of the macrocycle [1]. Formation of the in-plane complexes used to be enhanced by the addition of a small concentration of a metal ion with larger ionic radius (e.g. Cd2+, Hg2+, Pb2+) to the solution of the smaller one because the insertion of the larger metal ion into the ligand cavity is much faster. However, the OOP complex is considerably less stable. In its dome-distorted structure, two diagonal pyrrolic nitrogens are more accessible from the other side of the ligand, owing to the enhancement of their sp3 hybridization, hence, the metal center can be easily exchanged for the smaller one [1-2, 9]. This accessibility makes the realization of dinuclear out-of-plane monoporphyrins (2:1 complexes) possible if the metal ion possesses a low (single) positive charge and is large, i.e. its charge density is small enough, e.g. mercury(I), silver(I) and thallium(I) ions [1, 8-9]. Moreover, the out-of-plane position of the metal center, together with the dome-distorted structure (owing to the twisting of aromatic substituents from a nearly perpendicular position closer to the porphyrin plane) may promote the formation of so-called sandwich complexes of various compositions, in which two metal ions can coordinate to one macrocycle, and, reversely, one metal ion can concomitantly coordinate to two ligands, (Fig.3) [1-2, 9]. Lanthanide(III) ions form typical examples of metallo-oligoporphyrins because they are inclined to form complexes of higher coordination number (8-12). However, they are hard Lewis acids, hence, their insertion into the coordination cavity of the softer N-donor porphyrin ligand is a slow and complicated process in aqueous solutions. This Scheme 1. Simplified energy level diagram for the change of the porphyrin’s molecular orbitals in different types of complexes [9]. Figure 2. Absorption spectra of the free base (H2TSPP4–); the highly distorted, octabrominated free base (H2TSPPBr8 4–); a typical in-plane (PdIITSPP4–); and a typical out-of-plane metalloporphyrin (HgIITSPP4–) within the Soret range [1]. VALICSEK, KISS, FODOR, IMRAN, AND HORVÁTH Hungarian Journal of Industry and Chemistry 32 phenomenon partly originates from the stability of their aqua complexes. Due to the consequence of their Pearson-type hard character, they coordinate rather to the peripheral substituents of porphyrin (instead of the pyrrolic nitrogens), i.e. to the ionic group ensuring water-solubility if they possess similarly hard O-donor atoms (e.g. carboxy or sulfonatophenyl groups). At lower temperatures, under kinetic control, the early lanthanide(III) ions are not able to coordinate into the cavity, rather to the periphery; resulting in the formation of the tail-to-tail dimer of free-base ligands (as the tail used to be referred to as the periphery). Higher temperatures and thermodynamic control are also necessary for the insertion of metal ions into the cavity produced by four pyrrolic nitrogens; resulting in the formation of typical metalloporphyrin complexes. After the discovery of the possible coordination bonds between lanthanide ions and sulfonato substituents, the formation of lanthanide bisporphyrins may be realized as a tail-to-tail dimerization of two metallo- monoporphyrin complexes through a metal bridge; deviating from the head-to-head connection as in the case of typical sandwich complexes (head refers to the cavity; see Fig.3). On the basis of our previous experiences, the coordination position of lanthanide ions was influenced by the change in temperature [1-2], [11- 13]. During the investigation of the formation of “typical p-type hyperporphyrin” complexes (e.g. Tl+, Pb2+, Bi3+), the species possessing highly redshifted absorption bands are the end-products of metalation only in hydrophobic solvents, since they can appear in aqueous solutions as intermediates with shorter or longer lifetimes depending on the metal ion. The absorption spectra of the end-products of these transformation reactions are very similar (less redshifted) to those of the typical, common out-of-plane metallo-monoporphyrins (e.g. HgII-porphyrin in Fig.2). This phenomenon may be accounted for to the considerable coordination ability or the polarizing effect of water molecules, which can enable the complex to overcome the kinetic energy barrier towards the formation of the more stable structure, in which the metal center is located closer to the ligand plane, resulting in decreases in distortion, as well as redshift. Furthermore, “hyperporphyrins” can appear as intermediates in smaller amounts during the formation of typical, common out-of-plane metallo- monoporphyrins as well [1, 14]. In the case of “d-type hyperporphyrin” complexes (e.g. Mn3+, Co3+, Ni2+), the low-spin and ruffled complex with a contracted cavity can exist in a spin- isomerization equilibrium with the high-spin and planar forms, which not only exhibits less redshift, but rather blueshifted absorption bands compared to those of the free-base ligand. This reaction can be influenced by the strength of the M-N bonds (owing to the electronic effects of peripheral or axial substituents), due to the size of the coordination cavity (owing to the substitution or saturation of methylidene bridges or pyrrolic carbons) [1-2, 5-7]. 3.3. Photophysics Porphyrins represent one of the most interesting groups of compounds in terms of photophysical properties and biological significance. Due to their rigid structure and aromatic electronic system, they display two types of fluorescence: beside their relatively strong singlet-1 fluorescence in the range of 550–800 nm, weak and rare singlet-2 luminescence is observable between 400 and 550 nm upon excitation of the Soret band (Scheme 1) [1]. The quantum yields of S2-fluorescence are about 3 orders of magnitude lower than those of S1-fluorescence in the case of free-base porphyrins, especially ~1200- fold for H2TSPP4- (6.3×10-5 vs. 7.5%). However, this ratio decreases with metalation, mainly in the case of the formation of out-of-plane complexes. Since the structure of S2-excited porphyrins may be close to that of the dome-distorted OOP complexes that are already in the electronic ground state. Another consequence of this structural similarity (namely small Stokes shift) is that the directions of the shifts of S2-fluorescence bands invert (according to Soret absorption) between the in- plane (redshifted) and out-of-plane (blueshifted) complexes when compared to the free-base ligand [1]. Singlet-1 fluorescence bands exhibit blueshifts in both types (in-plane and out-of-plane) of metalloporphyrins, as a consequence of the aforementioned split in free-base ligands because of the presence of two protons, as well as their reduced symmetry (Scheme 1). Furthermore, almost all complexes exhibit similarly large Stokes shifts, as well as lifetimes and quantum yields. In the case of in-plane a) b) Figure 3. Potential structures of 3:2 bisporphyrin: (a) head-to-head or (b) tail-to-tail [2]. STUDY OF THE OUT-OF-PLANE METALLOPORPHYRINS 45(1) pp. 29–36 (2017) 33 metal centers, the spin-orbit coupling, as an electronic quenching effect, may be dominant. Whereas for typical out-of-plane metalloporphyrins, the distortion, as a steric effect, can enhance their non-radiative decay. The highly distorted (d- and p-type) hyperporphyrins, the paramagnetic in-plane complexes (e.g. FeIIITSPP3-), as well as the head-to-head-type OOP bisporphyrins {e.g. HgII 3(TSPP)2 6-} do not exhibit significant levels of luminescence at room temperature. Conversely, the paramagnetic out-of-plane complexes (e.g. LnIIITSPP3-) possess similar fluorescence properties to that of the diamagnetic ones because a paramagnetic metal ion can cause the disappearance of fluorescence by spin-orbit coupling only if it is located in the plane. In the OOP position, it is not able to perturb as efficiently the molecule orbitals of the macrocycle that result in the common absorption and emission out-of-plane characteristics [1-2, 9]. Lanthanide(III) bisporphyrins {LnIII 3(TSPP)2 3-} have many similarities in terms of absorption and emission properties to those of monoporphyrin complexes (LnIIITSPP3-). These may only originate from the very weak π–π interactions between the macrocycles in the tail-to-tail-type aggregations (Fig.3) [1-2, 11-13]. 3.4. Photochemistry Porphyrin derivatives are the main components of photosynthesis, synthetically as well. Since the overall quantum yield of fluorescence and intersystem crossing resulting in the formation of triplet states is in excess of 95%, merely a slight proportion of excitation energy is dissipated as heat from singlet states. This ratio is the major reason why porphyrins are efficient in terms of optical sensations and photosensitizations. Free-base and kinetically inert in-plane metalloporphyrins may be appropriate candidates to be applied in photocatalytic systems based on outer-sphere electron transfer. D-type hyperporphyrins can be particularly promising from this viewpoint owing to their distorted structure which may enhance the (photo)redox reactivity of these coordination compounds. In the presence of a suitable electron acceptor (methylviologen, MV2+) and donor (e.g. triethanolamine, TEOA), these complexes proved to be efficaciousl photocatalysts that transfer electrons between the ground-state reactants through an outer- sphere mechanism, generating the MV•+ radical cation. This system can be applied for the production of hydrogen from water [2, 5-7]. Contrarily, the inner-sphere photoredox reactions are characteristic of the out-of-plane metalloporphyrins because of this special coordination (Scheme 2): an irreversible photoinduced charge-transfer from the ligand to the metal center (ligand-to-metal charge transfer, LMCT) improves the efficiency of charge separation, which allows their utilization as catalysts in cyclic processes for the synthesis of chemicals capable of conserving light energy, hopefully in terms of the photochemical cleavage of water. Due to photoinduced LMCT the charge of the metal center decreases and its size increases, overall its charge density diminishes, hence, the coordinative bonds can easily split. The reduced metal ion can leave the cavity, primarily in polar solvents, and induce further redox reactions. The latter processes strongly depend on the stability of the reduced metalion in the actual medium. The oxidized and metal-free (cat)ionic radical of porphyrin is a very strong base: it is immediately protonated and forms the free-base radical, which is a long-lived and rather strong electron acceptor, especally in deaerated solutions. Since it would only oxidize water to oxygen at higher pHs, a slightly more efficent reducer (such as alcohols or aldehydes of low molecular weight) is needed, from which useful byproducts can be produced in terms of photocatalytic hydrogen generation. In the absence of any electron donor that promotes the regeneration of the porphyrin, it undergoes the primary photochemical processes; an overall four-electron oxidation involving a ring-cleavage, the end-product of which is a dioxo- tetrapyrrole derivative (bilindione). This ring-opening process can be followed by spectrophotometry owing to the disappearance of the Soret band, as well as the typical change in the region of Q bands [2, 4, 8-16]. Photochemical quantum yields of this ring-opening reaction (without regeneration) are about 2-3 orders of magnitude higher for the out-of-plane complexes (10-4 – 10-2) than for the free-base and in-plane metalloporphyrins (10-6 – 10-5). In addition, in the case of out-of-plane complexes, photoinduced dissociation in Scheme 2. Simplified demonstration of the mechanism for the inner-sphere photoredox reaction of an out-of- plane metalloporphyrin [8]. VALICSEK, KISS, FODOR, IMRAN, AND HORVÁTH Hungarian Journal of Industry and Chemistry 34 the absence of a redox reaction can occur, originating from their lability, and structural transformations to another complex form or conformer were observed in some cases as a photoinduced change of the type or measure of distortion (e.g. d- and p-type hyperporphyrins) [2, 4, 8-10, 14-16]. Besides the typical post-transition metal ions, lanthanide(III) ions were also applied for out-of-plane coordination because their contraction makes the fine- tuning of the out-of-plane distance possible, and their high negative redox potentials promote the photoinduced cleavage of water. Photochemical activities of their complexes confirm that the redox potentials of the metal centers are not the main determining factor, rather their out-of-plane distances [2, 11-13]. Deviating from the OOP complexes of post- transition metal ions, another stable photoproduct was observable during the photolysis of lanthanide(III) porphyrins. It displays a typical absorption band in the Q range (at ~600 nm), which may be assigned as a charge transfer between the metal ion and open-chain, dioxo-tetrapyrrole derivative (bilindione, see Scheme 2). Its oxo-groups, as donor atoms, may coordinate with the lanthanide ions, as a consequence of their similar Pearson-type hard characteristics, contrary to the softer post-transition metal ions [2, 11]. During the photolysis experiments, only small differences appeared between lanthanide(III) mono- and bisporphyrin complexes, which might confirm a special type of aggregation through the peripheral sulfonato substituents with weak π-π interactions (tail-to-tail, see Fig.3) [2, 11-13]. Deviating from these observations, the differences are much more significant in the case of the most typical, post-transition metallo-bisporphyrin compared to the monoporphyrin equivalent; namely between HgII 3(TSPP)2 6- and HgIITSPP4-. The overall quantum yield is ~2 orders of magnitude higher and the photoinduced dissociation of a metal ion became the dominant reaction in the head-to-head sandwich complex as a consequence of the strong π-π interactions [9-10]. 3.5. Quantum chemical calculations The main aims of our electronic structural calculations were to determine the primary consequences for the out- of-plane position of a metal center, and confirm the experimentally observed correlation between the UV- Vis spectral shifts and the coordination position of the metal center (in-plane or out-of-plane). In the light of these aspects, the unsubstituted porphin (H2P, C20H14N4) was used as a model for the calculations, instead of the tetrakis(sulfonatophenyl)porphyrin (H2TSPP4–, C44H26N4O12S4 4-). On the basis of the few comparative calculations that were conducted, the phenyl-, as well as the sulfonatophenyl substituents have negligible effects on the coordination of the metal center in the cavity. However, they can significantly influence the formation of bisporphyrin complexes, even in the case of head-to- head structures [4, 10]. According to our quantum chemical experience, the value of the critical radius became ~100 pm instead of the experimentally suggested ~75-90 pm as a consequence of the significant expansion of the coordination cavity to coplanarly incorporate the metal ions. The proportion of borderline cases, i.e. complexes with questionable structures (somewhere between in- plane and out-of-plane), increased with further post- transition metal ions (e.g. Ag2+ [15], Cd2+ [4], Tl3+ [16]) that possess ionic radii of ~90-95 pm. Calculated bond lengths (M-N) and atomic distances (N-N) considerably deviate from the expected ones supposed on the basis of the values of the deprotonated porphyrins (P2-). To describe this phenomenon, an axial ligand was applied to these metal centers to extract them out of the cavity. Consequently, expansion stopped, and the out-of-plane distance increased dramatically together with the degree of dome distortion and redshifts of absorption bands. From this point of view, two possible explanations can be supposed for the borderline-case complexes: the experimentally observed common OOP characteristics may originate from this expansion, tension; and small perturbations (e.g. the axial coordination in the calculation or photoexcitation in the experiments) may facilitate the metal center to adopt an out-of-plane position, too. Another possibility is that the method of calculation strongly prefers planar structures. In our time-dependent density functional theory (TD-DFT) calculations, the correlation found between the measured and calculated shifts associated with the position of the metal center was not totally linear, but nevertheless acceptable. The main exceptions were the borderline cases (high-spin Mn2+, Fe2+ and Zn2+) and the d-type hyperporphyrins (because their structures were determined to be totally planar), as well as the p-type hyperporphyrins (because a ruffled-like deformation did not superpose on their dome-like structure). The regression of correlation was much worse within the Soret band than in the case of the Q bands. The Soret band was also split in the calculations, which cannot be detected experimentally. On the basis of fur- ther experimental observations and doubts in the litera- ture, the validity of the theoretical model in use at pre- sent is questionable. Hence, the development of a more suitable one is in progress. 4. Conclusion In conclusion, it can be declared that the categorization of metalloporphyrins was complemented by the role of their distortion, which is primarily responsible for their spectral features, whereas the electronic structure of their metal centers is a secondary factor, with a considerable level of emphasis on the in-plane complexes. The position of the metal center (in-plane or out-of-plane) in the monoporphyrin complexes, as well as the type (head-to-head or tail-to-tail) of the bisporphyrin complexes can be determined on the basis of their UV-Vis absorption and emission properties. STUDY OF THE OUT-OF-PLANE METALLOPORPHYRINS 45(1) pp. 29–36 (2017) 35 Hyperporphyrin spectra can appear, owing to the peripheral substitution (octabromination) of free-base ligands. Furthermore, the high degree of redshift may disappear during the spin isomerization of d-type metalloporphyrins or the transformation of p-type ones. Consequently, the real origin cannot be an electronic but rather a steric effect, namely the measure of distortion, which can confirm the absence of their fluorescence. In terms of photochemical activity, several dissimilarities were found between the in-plane and out- of-plane metalloporphyrins; the most remarkable of them was the mechanism of their photoredox reactions: outer-sphere electron transfer is typical of the previous ones, while the inner-sphere equivalent is most prevalent for the latter ones. As a further consequence of the OOP position of the metal center, photoinduced dissociation and transformation reactions can occur within their complexes. In our electronic structural calculations, the number of borderline-case complexes expanded, on the basis of which common OOP characteristics that can be experimentally observed may acquire a novel explanation. Acknowledgement This research was supported by the Széchenyi 2020 Fund under the GINOP-2.3.2-15-2016-00016 and EFOP-3.6.1-16-2016-00015 projects. Assistance with quantum chemical calculations provided by Professor György Lendvay (Research Centre for Natural Sciences, Hungarian Academy of Sciences) is gratefully acknowledged. Finally, this manuscript is dedicated to the memory of Professor János Liszi, who, as the head of the Doctoral School for Chemistry at the University of Veszprém, praised the corresponding author’s PhD dissertation of a similar title in 2007. REFERENCES [1] Valicsek, Z.; Horváth, O.: Application of the elec- tronic spectra of porphyrins for analytical purpos- es: the effects of metal ions and structural distor- tions, Microchem. J., 2013 107, 47–62 DOI: 10.1016/j.microc.2012.07.002 [2] Horváth, O.; Valicsek, Z.; Fodor, M.A.; Major, M.M.; Imran, M.; Grampp, G.; Wankmüller, A.: Visible light-driven photophysics and photochem- istry of water-soluble metalloporphyrins, Coord. Chem. Rev., 2016 325, 59–66 DOI: 10.1016/j.ccr.2015.12.011 [3] Shelnutt, J.A.; Song, X.-Z.; Ma, J.-G.; Jia, S.-L.; Jentzen, W.; Medforth, C.J.: Nonplanar porphyrins and their significance in proteins, Chem. Soc. Rev., 1998 27, 31–41 DOI: 10.1039/A827031Z [4] Valicsek, Z.; Horváth, O.; Lendvay, G.; Kikaš, I.; Škorić, I.: Formation, photophysics, and photo- chemistry of cadmium(II) complexes with 5,10,15,20-tetrakis(4-sulfonatophenyl)porphyrin and its octabromo derivative: the effects of bro- mination and the axial hydroxo ligand, J. Photo- chem. Photobiol. A, 2011 218, 143–155 DOI: 10.1016/j.jphotochem.2010.12.014 [5] Fodor, M.A.; Horváth, O.; Fodor, L.; Grampp, G.; Wankmüller, A.: Photophysical and photocatalytic behavior of cobalt(III) 5,10,15,20-tetrakis(1- methylpyridinium-4-yl)porphyrin, Inorg. Chem. Commun., 2014 50, 110–112 DOI: 10.1016/j.inoche.2014.10.029 [6] Fodor, M.A.; Horváth, O.; Fodor, L.; Vazdar, K.; Grampp, G.; Wankmüller, A.: Photophysical and photochemical properties of manganese complexes with cationic porphyrin ligands: Effects of alkyl substituents and micellar environment, J. Photo- chem. Photobiol. A, 2016 328, 233–239 DOI: 10.1016/j.jphotochem.2016.06.011 [7] Major, M.M.; Horváth, O.; Fodor, M.A.; Fodor, L.; Valicsek, Z.; Grampp, G.; Wankmüller, A.: Photo- physical and photocatalytic behavior of nickel(II) 5,10,15,20-tetrakis(1-methylpyridinium-4-yl) por- phyrin, Inorg. Chem. Commun., 2016 73, 1–3 DOI: 10.1016/j.inoche.2016.09.001 [8] Horváth, O.; Valicsek, Z.; Harrach, G.; Lendvay, G.; Fodor, M.A.: Spectroscopic and photochemical properties of water-soluble metalloporphyrins of distorted structure, Coord. Chem. Rev., 2012 256, 1531–1545 DOI: 10.1016/j.ccr.2012.02.011 [9] Horváth, O.; Huszánk, R.; Valicsek, Z.; Lendvay, G.: Photophysics and photochemistry of kinetically labile, water-soluble porphyrin complexes, Coord. Chem. Rev., 2006 250, 1792–1803 DOI: 10.1016/j.ccr.2006.02.014 [10] Valicsek, Z.; Lendvay, G.; Horváth, O.: Equilibri- um, photophysical, photochemical and quantum chemical examination of anionic mercury(II) mono- and bisporphyrins, J. Phys. Chem. B, 2008 112(46), 14509–14524 DOI: 10.1021/jp804039s [11] Imran, M.; Szentgyörgyi, C.; Eller, G.; Valicsek, Z.; Horváth, O.: Peculiar photoinduced properties of water-soluble, early lanthanide(III) porphyrins, Inorg. Chem. Commun., 2015 52, 60–63 DOI: 10.1016/j.inoche.2014.12.016 [12] Kiss, M.P.; Imran, M.; Szentgyörgyi, C.; Valicsek, Z.; Horváth, O.: Peculiarities of the reactions be- tween early lanthanide(III) ions and an anionic porphyrin, Inorg. Chem. Commun., 2014 48, 22–25 DOI: 10.1016/j.inoche.2014.08.001 [13] Valicsek, Z.; Eller, G.; Horváth, O.: Equilibrium, photophysical and photochemical examination of anionic lanthanum(III) mono- and bisporphyrins: the effects of the out-of-plane structure, Dalton Trans., 2012 41, 13120–13131 DOI: 10.1039/C2DT31189E [14] Valicsek, Z.; Horváth, O.; Patonay, K.: Formation, photophysical and photochemical properties of wa- ter-soluble bismuth(III) porphyrins: the role of the charge and structure, J. Photochem. Photobiol. A, 2011 226, 23– 35 DOI: 10.1016/j.jphotochem.2011.10.011 VALICSEK, KISS, FODOR, IMRAN, AND HORVÁTH Hungarian Journal of Industry and Chemistry 36 [15] Harrach, G.; Valicsek, Z.; Horváth, O.: Water- soluble silver(II) and gold(III) porphyrins: the ef- fect of structural distortion on the photophysical and photochemical behavior, Inorg. Chem. Com- mun. 2011 14, 1756–1761 DOI: 10.1016/j.inoche.2011.08.003 [16] Valicsek, Z.; Horváth, O.: Formation, photophysics and photochemistry of thallium(III) 5,10,15,20- tetrakis(4-sulphonatophenyl)porphyrin: New sup- ports of typical sitting-atop features, J. Photochem. Photobiol. A, 2007 186, 1–7 DOI: 10.1016/j.jphotochem.2006.07.003 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 37–48 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0007 IMPROVEMENT POSSIBILITIES OF HETEROGENEOUS PHOTOCATALYSIS WITH THE AIM OF IN-FIELD USE ORSOLYA FÓNAGY, 1 PÉTER HEGEDŰS, 1 ERZSÉBET SZABÓ-BÁRDOS, 1 ANNAMÁRIA DOBRÁDI, 2 AND OTTÓ HORVÁTH 1* 1 Department of General and Inorganic Chemistry, University of Pannonia, Egyetem u. 10., H-8200 Veszprém, HUNGARY 2 Institute of Materials Engineering, University of Pannonia, Egyetem u. 10., H-8200 Veszprém, HUNGARY Heterogeneous photocatalysis can be successfully applied for the degradation of organic pollutants, although the efficiency of this method is insufficient. It can be increased by modification of the catalyst with precious metals (e.g. silver) or heterogeneous catalysis can be combined with other oxidative procedures such as ozonation. Another obstacle to the in- field use of the method is that separating the catalyst from the liquid phase is difficult. This problem can be eliminated by the immobilization of the catalyst. The poly(vinyl alcohol)-TiO2 immobilized catalyst has been prepared for this purpose. During TiO2-based photocatalysis, active oxygen species such as the hydroxyl radical, superoxide radical and hydrogen peroxide are produced. The formation rate of the oxidative •OH radicals was also determined in the case of the previously mentioned techniques. As a scavenger of this radical, coumarin was added. Keywords: TiO2-based photocatalysis; Triton X-100; Benzenesulfonic acid; Silver-deposited TiO2; Immobilization; Thermal decomposition of PVA 1. Introduction Nowadays, a great proportion of anthropogenic pollutants cannot be mineralized by traditional biological and physicochemical wastewater treatment procedures. Therefore, important developments in chemical wastewater treatment technologies and in their utilizations are required. For economic, technological and healthcare reasons, in terms of both type and quantity, it is important to use as few chemical additives as possible. Their applicability against remarkably different pollutants using minimal levels of energy consumption is even more important. Advanced Oxidation Processes satisfy these challenges [1]. A common characteristic of these processes is the generation of oxidative radicals, predominantly hydroxyl radicals, using solar radiation or other kinds of energy. Hydroxyl radicals are able to oxidize a great variety of organic compounds, therefore, heterogeneous photocatalysis has become an intensively studied field of research. In the pharmaceutical industry, benzenesulfonic acid is mostly used for producing other speciality chemicals. A variety of pharmaceutical drugs are prepared as salts of benzenesulfonic acid, known as besilates. Triton X-100 *Correspondence: horvath.otto@mk.uni-pannon.hu with an average n ≈ 9.5 is one of the most widely applied man-made non-ionic surfactants. Non-ionic surfactants are more stable than ionic tensides and not sensitive to pH and electrolytes. Besides the hydrophilic polyethoxy chain, it also contains a hydrophobic octylphenyl group. Both compounds can hardly be degraded by biological treatments [2–4]. Hence, a more efficient degradation method like heterogeneous photocatalysis is indispensable for the total mineralization of these pollutants. 1.1. Photodegradation of model compounds The titanium dioxide-mediated photocatalyzed degradation of benzenesulfonic acid (later abbreviated as BS) was investigated by monitoring the absorption and emission spectral changes, sulfate concentration, pH, as well as the total organic carbon (later abbreviated as TOC) content in both argon-saturated and aerated systems. During the degradation of benzenesulfonic acid, a characteristic change in the π→π* transition (typical of aromatic systems) in the absorption spectrum of BS (λmax=262 nm, ε=439 M–1 cm–1) can be observed. At first, there is a slight increase in the 260-270 nm range, while two new shoulders appear at around 280 and 300 nm and as the period of irradiation progresses the level of absorbance decreases. The reason for the increase in absorbance observed in the spectrum of this aromatic surfactant is that the absorbances of the intermediates FÓNAGY,HEGEDŰS,SZABÓ-BÁRDOS,DOBRÁDI, AND HORVÁTH Hungarian Journal of Industry and Chemistry 38 formed in the initial period of the process are higher than that of the starting compound. On the basis of the characteristic molar absorbances of 4- hydroxybenzenesulfonic acid (4-HBS) (λmax=277 nm, ε=868 M–1 cm–1) and 2,5-dihydroxybenzenesulfonic acid (λmax=302 nm, ε=3300 M–1 cm–1) as commercially available standards, the appearance of a shoulder at 280- 300 nm suggests the formation of hydroxylated intermediates. The results indicate that the initial step of the degradation is hydroxylation of the starting surfactant, resulting in the production of hydroxy- and dihydroxybenzenesulfonates [5–7]. Therefore, liquid chromatography-mass spectrometry analysis was utilized for the detection of intermediates. The sulfonate group of the substrate to be degraded can easily be oxidized by oxygen-containing reactants generated in these systems. Sulfate ions were detected in the irradiated reaction mixtures, which means that the hydroxylation reactions were accompanied by desulfonation. The initial rate of this process was lower in each case than that for the transformation (i.e. hydroxylation) of the starting compound (BS). These results indicate that the first step of mineralization is hydroxylation; desulfonation takes place afterwards. During the photocatalytic oxidation of BS, especially in the first 180 min, a continuous decrease of pH was observed in both the argon-saturated and the aerated systems, i.e. strong acidification occurred. Probably, there is a direct correlation between the formations of hydrogen and sulfate ions, due to the following reaction (Eq.(1)), which can take place between the sulfonate group and hydroxyl radical, the main oxidizing agent in these systems: •RSOHR•HSOOH•RSO 2 443   (1) The hydroxy species did not decay during the irradiation in the absence of dissolved oxygen. In the aerated system desulfonation and hydroxylation was much more efficient, moreover, a significant decrease in TOC took place during the initial stage. Further hydroxylation resulted in the cleavage of the aromatic system, through the formation of polyhydroxy derivatives, followed by ring-opening, which led to the production of aldehydes and carboxylic acids. Total mineralization was realized by the end of the photocatalysis. It has been proved that in this photocatalytic procedure the presence of dissolved oxygen is indispensable for the cleavage of the aromatic ring because hydroxyl radicals photochemically formed in the deaerated system alone are not able to break the C-C bond [5]. The degradation mechanism of Triton X-100 was also investigated [8]. In the case of this surfactant, hydroxyl radicals can attack in three ways; on the ethoxy chain, on the aromatic ring or on the alkyl chain. According to the absorption spectrum (no shift in the longer wavelength (275 nm) band and no new shoulders), no hydroxylation of the aromatic ring took place during the photocatalytic degradation. After the total disappearance of the starting surfactant (after approx. 60 min), an appreciable degree of absorbance arose at 275 nm, which indicates that aromatic intermediates were formed during the first hour. In liquid chromatographic experiments, Triton X-100 produced a multiple-peaked chromatogram. The retention time of these non-ionic surfactant components is in strong correlation with the length of the ethoxy chain. Intermediates formed from various lengths of ethoxy chains and others with aromatic rings were detected in the reaction mixture after irradiation started. The concentration of the shorter-chain molecules increased during the initial stage (in the first 10 min) of the irradiation, indicating that the fragmentation of the longer polyethoxy chains of the starting components initially increased the concentration of those with shorter ones, which was followed by a gradual decay. The component with a relative molecular weight of 206 Dalton could be clearly identified. Based on its structure, it was identified as the alkyl phenol part of the original Triton X-100. Its concentration peaked at around the 50th minute of irradiation. This is the result of the total cleavage of the polyethoxy chains without any oxidation of the rest of the original tenside molecules. Based on the results so far, it can be stated that hydroxyl radicals attack the ethoxy side-chain and subsequently the alkyl group, afterwards, the ring opens. 1.2. Challenges of heterogeneous photocatalysis The obstacle to the practical application of heterogeneous photocatalysis is that its efficiency is relatively low. Fortunately there are several ways to improve its mineralization effect: by the combination of heterogeneous catalysis with other oxidative procedures such as ozonation [6, 9-10] or by the modification of the catalyst with precious metals. The latter solution helps promote charge separation by accumulating the electrons on the surface of the silver particles [11-12]. The other obstacle is that the catalyst is in a suspended form, therefore, its separation from the liquid phase is difficult and makes the technology more expensive. This problem can be eliminated if the catalyst is immobilized. There are several methods for immobilization, for example, the catalyst particles can be fixed directly into polymers. Between a polymer and the catalyst usually a physical contact is formed, but in some cases, such as with poly(vinyl alcohol), a chemical bond can be formed. Unfortunately, this method is not free of disadvantages. The fact that immobilization decreases the active surface area of the catalyst has to be taken into consideration, thus, it will become less efficient [13–16]. Irradiation in the near UV range generates electron-hole pairs in the TiO2 nanoparticles. Under aerobic conditions the electrons in the conduction band are scavenged by adsorbed oxygen molecules with a relatively high quantum yield, producing O2• – ions, which are readily protonated in acidic media. Beides, IMPROVEMENT POSSIBILITIES OF HETEROGENEOUS PHOTOCATALYSIS 45(1) pp. 37–48 (2017) 39 the holes, in the absence of reducing species, oxidize the surface-adsorbed water or hydroxide ions to hydroxyl radicals. The •OH radical can also be formed by thermal reactions following the primary electron transfer processes that occur at the surface of semiconductor particles or in the liquid phase [13,17]. The hydroxyl radical is often assumed to be the major reactant responsible for the photooxidation of various organic molecules. The first step of the mineralization of benzenesulfonic acid is hydroxylation of the aromatic ring, while the degradation of Triton X- 100 is initialized by hydroxyl radicals, which attack the ethoxyl side chain (Fig.1). Due to the high reactivity and short lifetime of •OH, the direct detection of this species is difficult. Appropriate methods such as steady-state fluorescence were applied to measure the quantity of scavenged hydroxyl radicals. Recent studies have proved that several non- or weakly luminescent test molecules, such as coumarin, produce strongly luminescent compounds (7-hydroxycoumarin) with hydroxyl radicals. Thus, these molecules can be utilized for detecting and measuring •OH radicals produced by UV excitation of TiO2 particles in heterogeneous systems [18-19]. The main goal of our work was to investigate the effect of silver doping and immobilization of the catalyst applied for photocatalytic degradations. In this work, coumarin was used as a scavenger to determine the formation rate of hydroxyl radicals under various circumstances. 2. Experimental section 2.1. Materials For the experiments, the following materials were used: benzenesulfonic acid (C6H5O3S - Alfa Aesar), 4- hydroxybenzenesulfonic acid (C6H4(OH)O3S - Alfa Aesar), silver nitrate (AgNO3 - Aldrich), Degussa P25 TiO2: 25±5% rutile, 75±5% anatase with a specific surface area of 50 m2 g–1 (now called Evonik AEROXIDE® TiO2 P25), coumarin (C9H6O2 - VWR), Triton X-100 (C14H22O(C2H4O)9.5 - Alfa Aesar), and poly(vinyl alcohol) ((later abbreviated as PVA) M=146000-186000 g mol–1, hydrolysis: 99+% - Aldrich). The materials used for the HPLC and IC measurements were the following: methanol (CH3OH - Sigma-Aldrich), tetra-n-butylammonium bromide ((CH3CH2CH2CH2)4NBr - Sigma-Aldrich), sodium sulfate (Na2SO4 - VWR), 37% cc. hydrochloric acid (HCl - VWR) and sodium hydroxide (NaOH - VWR). High purity water, used in this study as a solvent, was double distilled and then purified using a Milli-Q system. 2.2. Analytics For the analysis, 4 cm3 samples were taken, the solid phase was removed by filtration using Millipore Millex- LCR PTFE 0.45 m filters. The absorption spectral changes of the reaction mixtures irradiated were followed using a PerkinElmer Lambda 25 UV-Vis spectrophotometer, while the emission spectra were recorded using a PerkinElmer LS 50B spectrofluorometer. Both measurements were carried out in 1 cm quartz cuvettes. Mineralization was followed by measuring the total organic carbon (TOC) concentration, utilizing a Thermo Electron Corporation TOC TN 1200 apparatus. The pH of the aqueous phase of the reaction mixture was determined by a SP10T electrode connected to a Consort C561 instrument. HPLC analyses of the benzenesulfonic acid samples were carried out using an Agilent 1100 instrument. A water:methanol (in a volume ratio of 95:5) solvent mixture containing 0.1 % (v/v) formic acid was the mobile phase with a flow rate of 1 cm3 min–1. During the experiments a 100×2 mm Synergi Hydro-RP C18 (Phenomenex, Torrance, CA, USA) column packed with 2.5 m particles was used. The column was thermostated at 30°C. The sample volume injected was 1 μl, photometric detection was applied at wavelengths of 262, 277 and 302 nm. Degradation of Triton X-100 was followed by the 1290 Infinity UHPLC system (Agilent Technologies Inc., Santa Clara, CA, USA). The mobile phase was 65:35 methanol:water for five minutes of analysis, then it was gradually increased to 75:25 over the following five minutes. The signal was detected with a DAD detector. The column used during the experiments was a 100×2 mm Synergi Hydro-RP C18 (Phenomenex, Torrance, CA, USA) column packed with 2.5 μm particles. The column was thermostated at 50 °C. The flow rate of the eluent was 1 cm3 min−1. During the study of the degradation mechanism an Agilent 6890N gas chromatograph with an Agilent DB- 5ms Ultra Inert column (30 m×0.25 mm×250 μm) was used for the separation. As a detector, an Agilent 5973 N-type mass spectrometer was used. The surface of immobilized catalysts was tested using a Philips/FEI XL30 scanning electron microscope at an accelerating voltage of 20 kV, using backscattered electron imaging. The composition was measured by an EDAX Genesis energy dispersive X-ray analyzer. Prior to testing, the samples were sputter coated by a 10 nm Figure 1. The role of •OH radicals in the degradation of the model compounds. FÓNAGY,HEGEDŰS,SZABÓ-BÁRDOS,DOBRÁDI, AND HORVÁTH Hungarian Journal of Industry and Chemistry 40 thick layer of gold. Elements of the qualitative composition do not overlap with the gold, therefore, the gold coating did not disturb the analysis. 2.3. Experimental circumstances The concentration of benzenesulfonic acid to be degraded was 10–3 M in each experiment, while that of the catalyst was 1 g dm–3 in all cases. For homogenization, the suspension was stirred for 20 minutes in the reactor before irradiation. The silver- deposited TiO2 was prepared by photoreduction from a mixture initially containing 10−4 M AgNO3 and 1 g dm–3 TiO2 as a photocatalyst. After 1 hour of photoreduction the catalyst, brown in color, was separated by vacuum filtration then dried at room temperature. To re-suspend the catalyst powder before each experiment an ultrasonic bath was applied in order to achieve a dispersion of particles sufficiently small in size. During subsequent ultrasonic treatment, the substrate to be degraded was added to the system only after a period of time, so that the •OH radicals generated during the sonication could recombine and not tamper with the experimental results. A detailed description of the production of immobilized catalyst can be found in our previous article [20]. The PVA-TiO2 catalyst was prepared by using a PVA with a molecular weight of 146 000-186 000 and a 99+% degree of hydrolysis and Degussa P-25 TiO2. Immobilization was carried out using a solution- casting method [21]. The PVA-TiO2 mixture was stirred at 90 then 60°C, subsequently a viscous solution was obtained, which was then poured into a petri dish and dried. The average weight of the immobilized catalyst was 0.47 g with a PVA content of 40%. Before its use, it was subjected to a 2 hour-long thermal treatment in an argon atmosphere to increase its stability. The Triton X- 100 concentration was 2∙10–4 M with 1.5 g dm−3 TiO2. It is important to note that when comparing the efficiencies of the catalyst in various forms, the amount of TiO2 immobilized in PVA (0.3 g) and used as a suspension in 200 cm3 of reaction mixture were both 1.5 g dm−3. Photochemical experiments were carried out using two different reactors and light-source setups. One was a laboratory-scale reactor with an effective volume of 2.5 dm3, equipped with a 40 W light tube located in the middle of the reactor. It emitted most of its energy at wavelengths exceeding 300 nm (λmax=350 nm, i.e. within the UVA range) [5]. The photon flux of the light source was I0 = 4.3∙10-6 mol photon dm–3 s–1, which was measured by the actinometry of trisoxalatoferrate(III). The heterogeneous reaction mixture was circulated by using a peristaltic pump through the reactor and a buffer vessel. Compressed air was bubbled through the reaction mixtures from gas bottles, to facilitate stirring and (with its O2 content) as an electron acceptor. O3 was produced by a LAB2B ozone generator, and introduced into the same airstream. The ozone concentration was determined by iodometry, using sodium iodide as a reagent and sodium thiosulfate for the titration of the iodine formed. The ozone dosage was estimated to be 0.35 mM min–1. The other setup consisted of one small reactor made of borosilicate glass (its volume was 250 cm3) [20]. The catalyst film was placed on a perforated glass platform (foil holder) standing on a glass frit, through which compressed air was bubbled from a gas bottle at a flow rate of 10 dm3 h−1, ensuring a constant oxygen concentration. The 200 cm3 reaction mixture (Triton X- 100 solution or distilled water in the case of pre- treatment) was circulated by a magnetic stirrer. It was irradiated from above using an Oriel LCS-100 solar simulator resulting in a light intensity of 72 mW cm–2 on the surface of the reaction mixture. Fig.2 shows the change in TOC during the photodegradation of benzenesulfonic acid under various experimental conditions. The band-gap energy of TiO2 is 3.2 eV, indicating that the catalyst can be excited by light of 360 nm in wavelength. The difference between the two light sources used is clearly visible. The radiation intensity of the solar simulator in the UV range is very low, so the total organic carbon content of the initial solution decreased by only 10 mg over 360 minutes, while total mineralization was achieved using the UV lamp within this period. 3. Results 3.1. Increase in photocatalytic efficiency using various methods As mentioned previously, the efficiency of heterogeneous photocatalysis can be increased in multiple ways: by combination with other oxidative procedures such as ozonation and by modification of the catalyst surface using noble metals. Since the first reaction in the mechanism of photocatalytic oxidative degradation of benzenesulfonic Figure 2. Change in the TOC concentration in the case of two different light sources used for irradiation. c(BS)0=10-3 M, air: 40 dm3 h-1, c(TiO2)=1 g dm-3 – air+TiO2+solar simulator  – air+TiO2+UV lamp IMPROVEMENT POSSIBILITIES OF HETEROGENEOUS PHOTOCATALYSIS 45(1) pp. 37–48 (2017) 41 acid is hydroxylation of the aromatic ring, the determination of the rates of formation of hydroxyl radicals as the main oxidant in these systems was indispensable as far as the elucidation of the mechanism under various circumstances was concerned. Coumarin was used as an •OH scavenger. An intensive band at 278 nm and a shoulder at around 300 nm is visible in the absorption spectrum of coumarin. The level of absorbance quickly decreased under irradiation (Fig.S1), while the concentration of the fluorescent 7-hydroxycoumarin rapidly increased during the initial period, and subsequently decreased gradually, i.e. degraded (Fig.S2). By depicting the intensity of the emission line at 453 nm as functions of the irradiation time, the initial rate of formation of 7- hydroxycoumarin, which is proportional to that of the hydroxyl radicals formed, can be determined from the initial gradients of the increasing sections, so it is easy to compare the ability of each process to generate hydroxyl radicals [18]. Before measuring the formation of 7- hydroxycoumarin in the photocatalytic oxidative systems that apply ozone, blind probes were used (Fig.3). The reaction between coumarin and ozone in the dark (at 10–4 M coumarin and 0.35 mM min-1 ozone concentrations) yielded a negligible spectral change after 2 hours. Although ozone itself is a strong oxidant, a very low rate of formation of 7-hydroxycoumarin (0.32 INT min–1) was observed under these experimental conditions. Upon UV irradiation of this probe, the rate of formation slightly increased (to 1.84 INT min–1). Ozone can be absorbed over a wide range, from infrared to vacuum UV, but is optimal at 254 nm. The light source applied in this work emits most of its radiation in the 300-380 nm range, where the molar absorbances of O3 are rather low. Thus, it can poorly promote the dissociation of ozone via the homolytic cleavage of a bond. Contrary to the blind probe with ozone, significant spectral changes were observed in irradiated reaction mixtures containing catalysts. The surface of the catalyst was also modified with silver and combined with ozonation, then the difference between the effects of the two methods of improvement was investigated by determining the rate of formation of oxidative •OH radicals (Fig.4). The combination of heterogeneous photocatalysis with ozonation resulted in a 2.5 fold increase in the rate of formation (from 7.42 INT min–1 to 19.4 INT min–1). A further increase was achieved by the application of surface-modified (silver-deposited) titanium dioxide (to 21.9 INT min–1). These results are in accordance with the fact that silver deposited on the surface of the catalyst diminishes the probability of the recombination of the photogenerated electron-hole pair, i.e. promotes charge separation [11]. 3.1.1. Degradation of benzenesulfonic acid Our research group has been studying the combination Figure 3. Change in the luminescence intensity of 7-hydroxycoumarin formed (λexc=332 nm). c(coumarin)0=10-4 mol dm-3, air: 40 dm3 h-1, c(O3)=0.35 mM, ℓ=1 cm  – O3  – O3+UV Figure 4. Change in the luminescence intensity of 7-hydroxycoumarin produced (λexc=332 nm) (a) and rates of formation of •OH radicals under various circumstances (b). c(coumarin)0=10–4 mol dm–3, air: 40 dm3 h–1, c(O3)=0.35 mM, c(catalyst)=1 g dm–3, ℓ=1 cm ♦ – air+TiO2+UV ● – O3+TiO2+UV ▲ – air+Ag-TiO2+UV a b FÓNAGY,HEGEDŰS,SZABÓ-BÁRDOS,DOBRÁDI, AND HORVÁTH Hungarian Journal of Industry and Chemistry 42 of heterogeneous photocatalysis and ozonation in detail. It has been established that a joint application of these procedures results in synergic effects [6]. Another possibility of enhancing the efficiency of photocatalysis is the modification of the catalyst surface by the deposition of silver. Our results in this work are compared with previous observations [17]. The spectral changes during the irradiations of suspensions containing Ag-TiO2 agree with earlier results. The degree of light absorption increases during the initial period (until the 50th minute), then it decreases (Fig.S3). The fine-band structure gradually disappears, while new bands arise at longer wavelengths (277 and 302 nm), suggesting the formation of isomers of hydroxy- and dihydroxybenzenesulfonic acid [6]. Characteristic changes can also be observed in the emission spectra of the samples in these experiments: an increase in the intensity of luminescence for 60 minutes, then a decrease accompanied by a red-shift of the 290 nm band (characteristic of the starting compound) and the appearance of a new emission band at about 330 nm (Fig.S4). The band-shift becomes considerable (46.5 nm) after the first hour of irradiation: from 286 nm to 333 nm. This finding also confirms the formation of hydroxy and dihydroxy intermediates with em.max.=305 nm and 352 nm, respectively. The changes in the concentrations of benzenesulfonic acid and the hydroxylated intermediates were followed by HPLC analyses (Fig.5). Of the three possible hydroxy isomers, the absolute concentration could only be determined for the commercially available 4-hydroxybenzenesulfonic acid (4-HBS, rt=0.86 min). Although the concentration of the substrate gradually decreases, the content of the benzenesulfonic acid in the reaction mixture is practically zero by the 180th minute of irradiation. On the basis of the concentration vs. time plot, the rate of decay of the model compound was determined by polynomial fit: 2.1∙10–2 mM min–1. Regarding the heterogeneous catalytic experiments (untreated TiO2 catalyst), 8∙10–3 mM min–1 was the degradation rate of the starting surfactant, while 1.4∙10–2 mM min–1 was measured in the case of combination with ozonation (Table 1) [6]. Mineralization of the model compound was followed by TOC (total organic carbon content) measurements. Taking the difference between the TOC value of the reaction mixture (TOC of the solution) and that of the benzenesulfonic acid (TOC of BS, calculated from its actual concentration), the TOC value regarding the intermediates in the solution can be obtained (Fig.6).This figure also displays the TOC values belonging to 4-HBS (calculated from its actual concentrations, see Fig.5) The concentration of the starting material and its TOC value rapidly decrease in the heterogeneous photocatalytic experiments involving silver-modified TiO2, while the amount of intermediates (predominantly hydroxy and dihydroxy derivatives) gradually increases. The initial rate of decrease in TOC for the reaction mixture was 1.58∙10–1 mg dm–3 min–1. The first step of the degradation of the model compound is hydroxylation. The highest initial rate of decay of BS was observed in the case of the system containing the silver-modified catalyst (Fig.7, Table 1). These results are in good agreement with those obtained for the experiments using coumarin, an efficient scavenger of •OH. The highest initial rate of formation of hydroxyl radicals was obtained on the Ag-TiO2 catalyst. Table 1. Comparison of the initial rates*. TOC Benzenesulfonic acid [mg dm-3 min-1] [mM min-1] air+ Ag-TiO2+UV 1.58∙10-1 2.3∙10-2 air+TiO2+UV 2.24∙10-1 6.5∙10-3 O3+TiO2+UV 5.92∙10-1 1.3∙10-2 *The error of the rate determinations is ±10%. Figure 5. Transformation of the initial compound and change in concentration of the 4-HBS intermediate. c(BS)0=10-3 M, air: 40 dm3 h-1, c(Ag-TiO2)=1 g dm-3 Figure 6. Change in TOC contents during the treatment. c(BS)0=10-3 M, air: 40 dm3 h-1, c(Ag-TiO2)=1 g dm-3  – TOC of the solution (measured)  – TOC of BS (calculated)  – TOC of intermediates (calculated)  – TOC of 4-HBS (calculated) IMPROVEMENT POSSIBILITIES OF HETEROGENEOUS PHOTOCATALYSIS 45(1) pp. 37–48 (2017) 43 The cleavage of the aromatic ring is indispensable in terms of the mineralization of the model compound, besides the fragmentation of the hydrocarbon chain. Szabó-Bárdos et al. rendered it probable that ring- opening takes place in the reaction of hydroxylated intermediates with reactive oxygen-containing species (O2• –/HO2•, O2( 1g)) [5]. Deviating from the decrease in the BS concentration, the lowest rate of the decrease in TOC, an indicator of mineralization, was observed in the case of the silver-modified catalyst, while the highest value (3.7 times higher than the previous one) was obtained for the combined procedure (Fig.7, Table 1). Based on our results, the application of heterogeneous photocatalysis with silver-modified TiO2 only accelerates the formation of hydroxyl radicals, while its combination with ozonation also increases the amount of other oxidative radicals (O2• –/HO2•, O2( 1g)). 3.2. Immobilization of the catalyst with poly(vinyl alcohol) The thermal treatment of the immobilized catalysts, the preparation of which was thoroughly described in Section 2.3, resulted in significant physical and chemical changes; on the one hand, their thickness varied (Fig.8), on the other hand, their color turned from their original white to brown (Fig.9). Song et al. described similar results [22]. The thickness of the foils prepared was tested by scanning electron microscopy prior to and after the thermal treatment. Following the treatment, the average thickness decreased from about 75.6 µm to about 69.6 µm (Fig.8). The shrinkage of about 6 µm can be attributed to the evaporation of water content. It has been established in our earlier studies that during the thermal treatment the polymer underwent dehydroxylation and cracking processes; compounds containing double and triple bonds as well as aromatic rings appeared on the surface of the foil, increasing its degree of light absorption [20]. Irradiation of the immobilized catalysts in the absence of the model compound led to the dissolution of the side-products of the cracking processes, and they also underwent heterogeneous photocatalytic degradation due to the presence of TiO2. Accordingly, stearic acid and its derivatives containing double bonds such as linoleic acid (9,12-octadecadienoic acid (Z,Z)- and oleic acid (9- octadecanoic acid) were detected by GC-MS analysis of Figure 7. Change in concentrations of BS (a) and TOC (b) under various circumstances. c(BS)0=10-3 M, air: 40 dm3 h-1, c(O3)=0.35 mM, c(catalyst)= 1 g dm-3 ♦ – air+TiO2+UV ● – O3+TiO2+UV ▲ – air+Ag-TiO2+UV a Figure 8. Cross sections of scanning electron micrographs of PVA-TiO2 foils before (a) and after (b) thermal treatment. a b b FÓNAGY,HEGEDŰS,SZABÓ-BÁRDOS,DOBRÁDI, AND HORVÁTH Hungarian Journal of Industry and Chemistry 44 the liquid phase. During the photocatalytic degradation, organic acids with shorter hydrocarbon chains were also produced, e.g. tetradecanoic acid (C14), caprylic acid (C8) and malic acid (C4) [20]. Although thermal treatment improves the stability of the foil, it restricts its efficiency because a considerable proportion of the oxidative radicals formed during photocatalysis are consumed in the reactions with the pollutants produced during the thermal treatment. Taking all these factors into consideration, how the temperature of the thermal treatment affects the stability of the immobilized catalyst was also studied. The foils heated at different temperatures were submersed in 200 cm3 of distilled water, then irradiated for 8 hours. The change in TOC of the liquid phase is shown in Fig.10/a, while the reductions in the foil masses are displayed in Fig.10/b. For the interpretation of the results, the increase in the temperature of the solution due to irradiation also has to be taken into account, because it enhances the solubility of poly(vinyl alcohol), i.e. the amount of the dissolved PVA, which increases the actual TOC value. The organic materials (PVA and the products of the thermal treatment) dissolved off the foil which also underwent photocatalytic degradation in the presence of TiO2. Both processes determine the change in TOC in the liquid phase, along with the mass reduction of the foil. The stability of the foils treated at lower temperatures (e.g. at 60oC) is rather modest, the PVA is dissolved in water, the TOC value of the liquid phase increased to 130 mg dm–3 by the end of the experiment, and the mass of the foil decreases considerably. The stability of the foils treated at higher temperatures (e.g. at 140oC or 200oC) significantly increased, while the TOC value of the liquid phase only grew slightly (to 34 mg dm–3 at 140°C, to 18 mg dm–3 at 200oC), in a similar fashion to the weight loss of the foils (140°C - m=93 mg, 200oC - m=24 mg). The results clearly indicate that by increasing the temperature of the thermal treatment, the stability of the foil is gradually improved, due to efficient crosslinking [21]. The thermally treated immobilized catalysts cannot be applied for photocatalytic purposes, because they undergo degradation themselves during the irradiation, increasing the TOC value of the solution phase. Besides, the brownish organic compounds formed in the cracking processes diminish the catalytic efficiency because they occupy the active sites on the catalyst surface. Hence, a compromise has to be made between stability and efficiency. Efforts were made to mineralize the products of cracking during a pre- treatment. Before the photocatalytic experiments, the thermally treated foils were cleaned by three sequential 8-hour-long pre-treatments, to remove the pollutants formed during the cracking processes. Since the pollutants formed whilst heated at the higher temperature (200oC) were rather difficult to remove, 140°C was applied for the thermal treatment of the foils in further experiments of ours [20]. 3.2.1. Efficiency of •OH production during pre- treatments of immobilized catalysts To detect the hydroxyl radicals formed during the irradiation of TiO2 particles embedded in the polymer, 10-4 M of coumarin as a scavenger was applied. Following the thermal treatment to enhance the stability of the foils, the undesirable products of this procedure have to be removed. In the first cycle of this process, when the foil is yet brown, the immobilized photocatalyst starts to work upon irradiation, but oxidative radicals formed in this process are mostly involved in the degradation of the products of cracking. The initial rate for the formation of 7-hydroxycoumarin just slightly exceeds the value observed in blind experiments (in the absence of a catalyst); 3∙10–2 INT min–1 (blind), 4∙10–2 INT min–1 (1st cycle) (Fig.11). Figure 9. Immobilized catalysts before and after heat treatment. Figure 10. Change in the TOC content of the liquid phase after various heat treatments (a) and the weight loss of foils after 8 hours of irradiation (b) 200 cm3 distilled water, air: 10 dm3 h-1,-c(TiO2)=1.5 g dm-3 IMPROVEMENT POSSIBILITIES OF HETEROGENEOUS PHOTOCATALYSIS 45(1) pp. 37–48 (2017) 45 The rate of formation of the hydroxyl radicals significantly increases during the 2nd cycle; the emission intensity at the end of this cycle (v0=10–1 INT min–1) is more than twice as high as in the previous one. The results obtained during the 3rd cycle are similar to those during the 2nd one. Practically, during the three sequential 8-hour-long periods of irradiation, the polluting compounds that originate from the thermal treatment are totally degraded, hence, the active sites on the surface of the catalyst become more accessible, resulting in a much higher rate for the formation of oxidative radicals (in the 4th cycle v0=2.5-2.8∙10–1 INT min–1, see Fig.11). In other applications, the radical producing efficiency of the catalyst does not change significantly, i.e. the foil can be applied at constant efficacy. The lower value determined during the 7th cycle suggests a slow rate of degradation of the immobilized catalyst, which has to be taken into consideration during its application. 3.2.2. Degradation of the model compound The photocatalytic degradation of the model compound was investigated following the cleaning procedure. The average degree of degradation of Triton X-100 in three sequential cycles was 62%; the deviations might originate from measurement errors (Fig.12). This value is significantly lower than that observed in the case of the suspended catalyst (90%). Although the catalyst concentrations with regard to the volume of solution were equal (1.5 g dm–3), immobilization considerably diminished the specific surface of the catalyst, leading to a substantial decrease in its efficiency. Mineralization of Triton X-100 takes place via cleavage-producing intermediates. The cleavage of the alkyl-phenyl part of the starting components of Triton X-100 took place from the beginning of the irradiation period. As completive intermediates originated from such a cleavage, ethylene glycol and dioxolane derivatives were detected by GC-MS. At the beginning of the irradiation the edge of the starting molecules or the initially generated intermediates mineralize to CO2. Therefore, the decrease in TOC is much less than expected from the concentration change. The concentration changes of Triton X-100 components and the typical intermediates are in accordance with our previous results using a similar system that applied a suspended catalyst (see section 1.1). Remarkably, GC- MS analyses did not detect any intermediates which might derive from the PVA support. This suggested that three pre-treatment cycles of irradiation and rinsing were sufficient for a satisfactory level of removal of the derivatives formed during the thermal treatment. 3.2.3. Surface morphology of foils - SEM analysis The behavior of the active surface of immobilized photocatalysts after repeated use raises an important issue. The surfaces of titanium dioxide containing PVA before and after thermal treatment are shown in Figs.13/a and 13/b, respectively. Some TiO2-rich islands (Fig.14) are not fully covered by the PVA film, these are visible on the surface. Otherwise, no significant changes are observed after the thermal treatment. However, simulated solar radiation (repeated usage) induces visible changes in the morphology of the b Figure 11. Change in the intensity of luminescence of 7-hydroxycoumarin formed. c(coumarin)0=10-4 M, air: 10 dm3 h-1, λexc=332 nm Figure 12. Change in concentration change of Triton X-100 during the irradiation. c(TX-100)0=2∙10-4 M, air: 10 dm3 h-1, c(TiO2)=1.5 g dm-3 FÓNAGY,HEGEDŰS,SZABÓ-BÁRDOS,DOBRÁDI, AND HORVÁTH Hungarian Journal of Industry and Chemistry 46 surface. Such a period of radiation partially or totally removes the PVA layer covering the TiO2 agglomerates. This process is enhanced after prolonged use: after four cycles the homogeneously distributed particles are already visible (Fig.13/c). Extended use - seven cycles - will enhance grain boundaries by removing more of the PVA matrix, indeed some of the particles become loosely bounded to the rest of the agglomerate (Fig.13/d). The observed changes in thickness strongly correlate with the degradation of films (Fig.15). After four cycles, a significant decrease in thickness can be observed: the foil thickness decreases to about 58.5 µm (from 69.6 µm). This corresponds to a decrease of about 2.7 µm / cycle. As for the thickness measured after seven cycles, the last three cycles only represent a decrease of about 0.9 µm/cycle. These observations imply a more intense degree of degradation during the initial stage. 4. Conclusion Heterogeneous photocatalytic procedures are widely applied for the degradation of organic pollutants. The degradation of benzenesulfonic acid and Triton X-100 by heterogeneous photocatalysis was investigated under various circumstances. The formation of hydroxyl radicals as one of the major reactants responsible for the photooxidation of these molecules was also measured in the case of both setups of the catalyst. Silver deposition on the surface of the TiO2 catalyst significantly increased the initial degradation (hydroxylation) rate of BS, due to the hindered recombination of the photogenerated electron-hole pairs, increasing the rate of •OH production. However, mineralization of the intermediates was not accelerated by this method, deviating from the combination of the TiO2-based photocatalysis with ozonation. The degradation efficiency in the latter procedure indicated a synergic effect on the mineralization rate as a consequence of the enhanced production of oxidative agents other than •OH radicals [6]. These reactive oxidants play key roles in the ring-opening processes needed for the mineralization of aromatic intermediates. Quantitative determinations of such reactants (e.g. O2• – and O2( 1g)) to gain a deeper understanding of the mineralization mechanism are in progress. Our results indicate that the methods applied for the enhancement of the degradation efficiency in heterogeneous photocatalysis affect the generation of the oxidizing agents in different ways, which ought to be taken into consideration in further studies. In general, suspensions of the photocatalyst are used in these procedures. Thus, at the end of these processes the semiconductor particles have to be separated from the cleaned solution. This operation can be taken out by immobilization of the catalyst, which was achieved by using poly(vinyl alcohol) in this work. Figure 13. Scanning electron micrographs of foil surfaces. a) before and b) after thermal treatment c) after 4 irradiation cycles d) after 7 cycles of irradiation a b c d Figure 14. Energy dispersive X-Ray spectra of the TiO2 island. In te ns it y / au IMPROVEMENT POSSIBILITIES OF HETEROGENEOUS PHOTOCATALYSIS 45(1) pp. 37–48 (2017) 47 The stability of the PVA support proved to be too low for practical use. The stability of the PVA-based foils could be considerably improved by thermal treatment. However, before photocatalytic application the water- soluble products of the heating process had to be removed from the surface of the PVA-TiO2 foil. According to the efficiency of •OH formation, three cycles of irradiation and rinsing proved to be sufficient for this purpose. The foil pre-treated by this method could be applied for the photocatalytic degradation of a frequently used non-ionic detergent, Triton X-100. When using the immobilized catalyst, it should be noted that the foil catalyst is degraded, its thickness decreases and the catalyst particles become more and more accessible on the surface. Acknowledgements This work was supported by Széchenyi 2020 under the projects GINOP-2.3.2-15-2016-00016 and EFOP-3.6.1- 16-2016-00015. REFERENCES [1] Dewil, R.; Mantzavinos, D.; Poulios, I.; Rodrigo, M. A.: New perspectives for Advanced Oxidation Processes, J. Environ. Manage., 2017 195 93–99 DOI: 10.1016/j.jenvman.2017.04.010 [2] Hack, N.; Reinwand, C.; Abbt-Braun, G.; Horn, H.; Frimmel, F. H.: Biodegradation of phenol, salicylic acid, benzenesulfonic acid, and iomeprol by Pseudomonas fluorescens in the capillary fringe, J. Contam. Hydrol., 2015 183 40–54 DOI: 10.1016/j.jconhyd.2015.10.005 [3] Mohan, P. K.; Nakhla, G.; Yanful, E. K.: Biokinetics of biodegradation of surfactants under aerobic, anoxic and anaerobic conditions, Water Res., 2006 40(3) 533–540 DOI: 10.1016/j.watres.2005.11.030 [4] Wyrwas, B.; Dymaczewski, Z.; Zgoła- Grześkowiak, A.; Szymański, A.; Frańska, M.; Kruszelnicka, I.; Ginter-Kramarczyk, D.; Cyplik, P.; Ławniczak, Ł.; Chrzanowski, Ł.: Biodegradation of Triton X-100 and its primary metabolites by a bacterial community isolated from activated sludge, J. Environ. Manage., 2013 128 292–299 DOI: 10.1016/j.jenvman.2013.05.028 [5] Szabó-Bárdos, E.; Markovics, O.; Horváth, O.; Törő, N.; Kiss, G.: Photocatalytic degradation of benzenesulfonate on colloidal titanium dioxide, Water Res., 2011 45(4) 1617–1628 DOI: 10.1016/j.watres.2010.11.045 [6] Zsilák, Z.; Szabó-Bárdos, E.; Fónagy, O.; Horváth, O.; Horváth, K.; Hajós, P.: Degradation of benzenesulfonate by heterogeneous photocatalysis combined with ozonation, Catal. Today, 2014 230 55–60 DOI: 10.1016/j.cattod.2013.10.039 [7] Zsilák, Z.; Fónagy, O.; Szabó-Bárdos, E.; Horváth, O.; Horváth, K.; Hajós, P.: Degradation of industrial surfactants by photocatalysis combined with ozonation, Environ. Sci. Pollut. Res., 2014 21(19) 11126–11134 DOI: 10.1007/s11356-014-2527-2 [8] Hegedűs, P.; Szabó-Bárdos, E.; Horváth, O.; Horváth, K.; Hajós, P.: TiO2-mediated photocatalytic mineralization of a non-ionic detergent: Comparison and combination with other advanced oxidation procedures, Materials (Basel)., 2015 8(1) 231–250 DOI: 10.3390/ma8010231 [9] Rodríguez, E. M.; Márquez, G.; León, E. A.; Álvarez, P. M.; Amat, A. M.; Beltrán, F. J.: Mechanism considerations for photocatalytic oxidation, ozonation and photocatalytic ozonation of some pharmaceutical compounds in water, J. Environ. Manage., 2013 127 114–124 DOI: 10.1016/j.jenvman.2013.04.024 [10] Agustina, T. E.; Ang, H. M.; Vareek, V. K.: A review of synergistic effect of photocatalysis and ozonation on wastewater treatment, J. Photochem. Photobiol. C Photochem. Rev., 2005 6(4) 264–273 DOI: 10.1016/j.jphotochemrev.2005.12.003 Figure 15. Scanning electron micrographs of cross sections of film. a) after 4 cycles of irradiation b) after 7 cycles of irradiation. a b FÓNAGY,HEGEDŰS,SZABÓ-BÁRDOS,DOBRÁDI, AND HORVÁTH Hungarian Journal of Industry and Chemistry 48 [11] Szabó-Bárdos, E.; Czili, H.; Horváth, A.: Photocatalytic oxidation of oxalic acid enhanced by silver deposition on a TiO 2 surface, J. Photochem. Photobiol. A., 2003 154 195–201 DOI: 10.1016/S1010-6030(02)00330-1 [12] Ma, J.; Guo, X.; Zhang, Y.; Ge, H.: Catalytic performance of TiO2@Ag composites prepared by modified photodeposition method, Chem. Eng. J., 2014 258 247–253 DOI: 10.1016/j.cej.2014.06.120 [13] Chong, M. N.; Jin, B.; Chow, C. W. K.; Saint, C.: Recent developments in photocatalytic water treatment technology: A review, Water Res., 2010 44(10) 2997–3027 DOI: 10.1016/j.watres.2010.02.039 [14] Gar Alalm, M.; Tawfik, A.; Ookawara, S.: Enhancement of photocatalytic activity of TiO2 by immobilization on activated carbon for degradation of pharmaceuticals, J. Environ. Chem. Eng., 2016 4(2) 1929–1937 DOI: 10.1016/j.jece.2016.03.023 [15] Bouarioua, A.; Zerdaoui, M.: Photocatalytic activities of TiO2 layers immobilized on glass substrates by dip-coating technique toward the decolorization of methyl orange as a model organic pollutant, J. Environ. Chem. Eng., 2017 5(2) 1565– 1574 DOI: 10.1016/j.jece.2017.02.025 [16] Veréb, G.; Ambrus, Z.; Pap, Z.; Mogyorósi, K.; Dombi, A.; Hernádi, K.: Immobilization of crystallized photocatalysts on ceramic paper by titanium(IV) ethoxide and photocatalytic decomposition of phenol, React. Kinet. Mech. Catal., 2014 113(1) 293–303 DOI: 10.1007/s11144-014- 0734-y [17] Kondrakov, A. O.; Ignatev, A. N.; Lunin, V. V.; Frimmel, F. H.; Bräse, S.; Horn, H.: Roles of water and dissolved oxygen in photocatalytic generation of free OH radicals in aqueous TiO2 suspensions: An isotope labeling study, Appl. Catal. B Environ., 2016 182 424–430 DOI:10.1016/j.apcatb.2015.09.038 [18] Czili, H.; Horváth, A.: Applicability of coumarin for detecting and measuring hydroxyl radicals generated by photoexcitation of TiO2 nanoparticles, Appl. Catal. B Environ., 2008 81(3–4) 295–302 DOI:10.1016/j.apcatb.2008.01.001 [19] Ishibashi, K. I.; Fujishima, A.; Watanabe, T.; Hashimoto, K.: Detection of active oxidative species in TiO2 photocatalysis using the fluorescence technique, Electrochem. commun., 2000 2(3) 207–210 DOI:10.1016/S1388-2481(00)00006-0 [20] Hegedűs, P.; Szabó-Bárdos, E.; Horváth, O.; Szabó, P.; Horváth, K.: Investigation of a TiO2 photocatalyst immobilized with poly(vinyl alcohol), Catal. Today., 2017 284 179–186 DOI:10.1016/j.cattod.2016.11.050 [21] Lei, P.; Wang, F.; Gao, X.; Ding, Y.; Zhang, S.; Zhao, J.; Liu, S.; Yang, M.: Immobilization of TiO2 nanoparticles in polymeric substrates by chemical bonding for multi-cycle photodegradation of organic pollutants, J. Hazard. Mater., 2012 227–228 185–194 DOI:10.1016/j.jhazmat.2012.05.029 [22] Song, Y.; Zhang, J.; Yang, H.; Xu, S.; Jiang, L.; Dan, Y.: Preparation and visible light-induced photo-catalytic activity of H-PVA/TiO2 composite loaded on glass via sol-gel method, Appl. Surf. Sci., 2014 292 978–985 DOI: 10.1016/j.apsusc.2013.12.090 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 49–59 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0008 ADAPTING THE SDEWES INDEX TO TWO HUNGARIAN CITIES VIKTOR SEBESTYÉN, VIOLA SOMOGYI,* SZANDRA SZŐKE, AND ANETT UTASI Institute of Environmental Engineering, University of Pannonia, 10 Egyetem str., Veszprém, H-8200, HUNGARY Numerous cities aim to mitigate their contribution to climate change and provide a liveable environment in the context of sustainable development. In order to measure these efforts, benchmarking performance would be a good solution. Methods for environmental analysis have their limitations when it comes to evaluating a city and other aggregated indicators focus on certain aspects of a sustainable or liveable settlement. The SDEWES In- dex was used for benchmarking several cities of different sizes in terms of metrics related to energy, water and environmental systems successfully thus it was chosen to compare the performance of Veszprém and Zalae- gerszeg, two environmentally conscious Hungarian county seats of roughly the same size and population. The SDEWES Index consists of 7 dimensions, namely energy consumption, industrial profile with CO2 emissions, CO2-saving measures, R&D, renewable energy potential and utilization, water and environmental quality, and social environment and sustainability policy. Each dimension is composed of 5 indicators that provide infor- mation on sustainable development of energy, water and environmental systems in cities. Using the SDEWES Index the strengths and weaknesses of the two cities are highlighted, locating those key parameters where im- provement can be achieved. Both for Veszprém and Zalaegerszeg progress could be realized concerning ener- gy-saving measures and the proportion of green areas could be increased. To improve the method and facilitate a more comprehensive comparison of cities of differing sizes, data should be provided concerning the territory or population. Also, the definition and inclusion of a worst and best case scenario that takes into account the pa- rameters would be advantageous in terms of a comparison. These were named ‘horror’ and SDEWES cities by the authors, respectively. Keywords: SDEWES Index, Sustainability, City development, Sustainable Energy Action Plan, City sample 1. Introduction Sustainability is a key issue when it comes to the development of cities. In 2014, 54% of the world’s population lived in urban areas and according to the prognosis, the proportion will be as high as 66% by 2050 [1]. While the number and population of megacities is on the rise, 43% of urban dwellers lived in settlements consisting of less than 300,000 inhabitants in 2014 (in Europe the corresponding data was 58%) and only a modest decrease is estimated by 2030 [2]. Looking at these numbers, it is easy to notice that urbanised areas have a huge impact on achieving sustainable development. Various cities have started to address this issue and several indicators and comparisons were created to measure specific aspects of sustainable development. A few of these numerous examples are listed in this paper. The City Development Index [3] studies the municipalities from social and governance aspects. The Global Power City Index (GPCI) has ranked 40 metropolises since 2008 [4] considering 70 individual indicators regarding the environment, liveability and *Correspondence: somogyiv@uni-pannon.hu R&D among others. The Green City Index is focused on the environmental sustainability of large cities [5]. Carbon footprints of twelve metropolitan areas [6] and the San Francisco Bay Area [7] were determined. Several others are listed by López-Ruiz et al. [8] but usually smaller towns do not fall within the scope of these benchmarks. Environmental analysis, on the other hand, facilitates the evaluation of the impacts of human activities (different actions, projects or investments) with regard to the local environment, economy and society. In this way it provides information on the status quo and helps the practical implementation of sustainable development by focusing attention on the points to be improved [9]. Several methods have been developed to carry out the procedure: checklists [10], the matrix technique [11], the network approach [12], GIS-based methods [13] and quantitative methods [14] may be used to evaluate environmental impacts. Aggregating methods such as the Global Pollution Index (IPG) [15] may be suitable up to a point in providing a comprehensive sustainability analysis as they are only based on immission values. Several multiple-criteria decision-making (MCDM) techniques designed to assist with decision-making, e.g. the Analytic Hierarchy Process (AHP) [16] or the Technique for Order Preference by Similarity to Ideal Solution systems (TOPSIS) [17] combined with Simple SEBESTYÉN, SOMOGYI, SZŐKE, AND UTASI Hungarian Journal of Industry and Chemistry 50 Additive Weighting (SAW) [18], can be used as well [19]. These hierarchical methods (TOPSIS, SAW and AHP) rank the examined parameters which may be useful when deciding between options of individual investments but are problematic in terms of adapting them to the decision-making process with regard to development strategies of the settlements. The Sustainable Development of Energy, Water and Environment Systems (SDEWES) City Sustainability Index was developed to overcome the disadvantages and limitations of other measures with regard to benchmarking the performance of cities in terms of energy, water and environment systems. So far a list of 58 cities assessed by the SDEWES Index can be accessed online on the SDEWES Centre homepage [20]. Also, articles concerning the benchmarking of 12 South East European cities (such as Athens and Belgrade) [21], 22 Mediterranean port cities (e.g. Barcelona and Venice) [22], and a further 18 South East European cities (including Budapest and Pécs) [23] were published, and the inventory will no doubt be expanded upon in the near future. In this paper two cities were evaluated by using the composite SDEWES Index and to test the method itself. Veszprém and Zalaegerszeg are two Hungarian county seats of roughly the same size and population, both are environmentally conscious, and are aiming to become environmentally friendly, liveable and sustainable cities. Veszprém is near Lake Balaton with a population of around sixty-two thousand people and and a surface area of 126.9 km2. The city has won the Climate Star award † [24] and aims to become an eco-city. In its Energy Strategy [25] the following objectives were set by 2026:  20% of the energy demand should be satisfied by renewable energy resources while the energy renovation of public and residential buildings should be 70% complete resulting in a reduction in greenhouse gas (GHG) emissions of 25%;  35% of transportation has to be conducted by means of public transport with environmentally friendly vehicles that are less than 10 years old and 10% of the vehicle-kilometres should be undertaken by bicycles;  increasing the proportion of green areas to 25 m2/capita and 60% of rainwater should be reused in some way. Zalaegerszeg is situated in the west of the country and consists of sixty-two thousand inhabitants and a surface area of 100 km2. Since the millennium its urban development strategies have focused on becoming a “Sustainable City”. In the strategy formulated in 2014 † The Climate Star award was founded by the Climate Alliance with the aim of demonstrating how climate protection initiatives can be implemented from the grass roots up [24]. Cities with initiatives in the fields of sustainable energy, mobility, consumption, urban and regional development, and citizen involvement may apply for the call in four categories. [26] a major goal was to improve energy efficiency by 20% while producing more than 20% of its energy using local renewable resources by 2030. This would result in a reduction of 36% in terms of energy-related costs. Taking 2005 as a base the GHG emissions should be reduced by 20% by 2030 while the particle pollution PM10 is planned to be mitigated by 10% by 2023. Besides achieving these indicators, the city council aims to create and strengthen its image of being an environmentally conscious, modern and sustainable green city. 2. Experimental The SDEWES Index consists of 7 dimensions and 35 main indicators (Table 1). The indicators of each dimension are explained in detail in Ref. [21]; only those that need further clarification or some sort of adjustment due to problems concerning the access of data are highlighted in this paper. The data for each indicator were normalised according to the Min-Max method [27]. Depending on whether the lower or higher values are more desirable, either Eq. (1) or (2) is used [21-22]. An example of the first case, i.e. when lower values are favourable, would be CO2 emissions, while the normalised data for the number of local universities would be calculated by Eq. (2). Since the leader (��,����� = max���,��) equals 1 and the laggard (��,����� = min���,��) 0, if only two cities are compared and the values are identical, the denominator would become 0. To avoid this ��,����� should be set to 0 for such cases. Another solution would be to include further cities in the benchmark. ��,����� = ���,������max���,��� �min���,���max���,��� (1) ��,����� = ���,������min���,��� �max���,���min���,��� (2) where: I – normalised value of the indicator, x – dimension number, y – indicator number within a dimension, Cj – jth city, i – data input before normalization. Value aggregation is done according to SDEWES���� = ∑ ∑ �� � ��� ��,����� � ��� , (3) where ∑ �� = 1� ��� and αx is the weight of the xth dimension. The SDEWES Index of the jth city is calculated by a double summation, where α1 and α5 are 0.22 since these dimensions include energy and CO2 emissions data. Other dimensions are weighed less (αx=0.11) as they are not directly related to the sustainable energy action plan [21-22]. ADAPTING THE SDEWES INDEX 45(1) pp. 49–59 (2017) 51 Table 1. The Dimensions and Indicators [21-23]. Dimensions D1: Energy Consumption and Climate D2: Penetration of Energy and CO2-Saving Measures D3: Renewable Energy Potential and Utilization D4: Water and Environmental Quality D5: CO2 Emissions and Industrial Profile D6: City Planning and Social Welfare D7: R&D, Innovation and Sustainability Policy In di ca to rs Energy consumption of buildings [MWh] Sustainable Energy Action Plan (SEAP) Solar energy potential [Wh/m2/day] Domestic water consumption [m3/capita] CO2 emission of buildings [t CO2] Price of a public transport ticket [EUR] R&D and innovation policy orientation Energy consumption of transport [MWh] Combined heat and power-based district (H/C) Wind energy potential [m/s] Water quality index [/100] CO2 emissions of transport [t CO2] Urban form and protected sites National patents in clean technologies Total energy consumption per capita [MWh/capita] Energy savings in end-usage (buildings) Geothermal energy potential [mW/m2] Annual mean PM10 concentration [µg/m3] Average CO2 intensity [t CO2/MWh] GDP per capita [PPP$ national] Local public/private universities Heating Degree Days (HDD) [day °C] Density of the public transport network Renewable energy usage for electricity [%] Ecological footprint [gha/capita] Number of CO2-intense industries Inequality adjusted well- being (HPI) National h- index of scientific publications Cooling Degree Days (CDD) [day °C] Efficient public- lighting armatures Biofuel share in transport [%] Biocapacity [gha/capita] Airport Carbon Accreditation (levels) Tertiary education rate (national) Reduction Target for CO2 emission reduction (2020) 2.1. Data of Veszprém and Zalaegerszeg As of October 2017, neither of the cities are signatories of the Covenant of Mayors movement, therefore, alternative sources of data had to be found. Besides the sources suggested by the developer of the index, the energy and integrated city development strategies were used to retrieve data. Necessary changes and the simplification of the original method is explained in detail below. 2.1.1. Energy Consumption and Climate (D1): The energy consumption of buildings (municipal, residential and commercial) and transportation (public, private and the municipal vehicle fleet) are indicators on their own (Table 2) but also included in terms of the total energy use per capita [21]: � �(��) = �∑ �b � ��� �∑ �� � ��� ��g��d� ����� (4) where  E – total energy consumption (MWh),  P(Cj) – population of the jth city (capita),  Eb – energy consumption of buildings (1: municipal, 2: residential and 3: commercial) (MWh),  Et – energy consumption of transport (1: public, 2: private and 3: municipal vehicle fleet) (MWh),  Eg – energy consumption of public lighting (MWh),  Ed – energy consumption of industry (MWh). The energy consumption of transport was calculated based on the number of vehicles registered according to the energy strategies of the cities [25-26] by presuming an average mileage of 15,000 km/year and average consumption of 7.5 l/100 km [25]. The energy content of diesel and gasoline was assumed to be 10.83 kWh/l and 8.89 kWh/l, respectively. Data for commercial buildings are only included in the total energy consumption indicator as the consumption of the service sector and industries was not collected separately by the cities. As for the municipal vehicle fleet, due to a lack of data for Veszprém, this had to be neglected for both towns. It has to be noted that the data for Veszprém were from between 2007 and 2009 as stated in the strategy while for Zalaegerszeg information SEBESTYÉN, SOMOGYI, SZŐKE, AND UTASI Hungarian Journal of Industry and Chemistry 52 could only be obtained from between 2012 and 2013. In later documents, only improvements are mentioned, newly obtained data on overall consumption is not stated. Monitoring the energy consumption of public buildings and lighting is an issue for both cities that needs to be solved. 2.1.2. Penetration of Energy and CO2-saving Measures (D2): Neither of the cities have a Sustainable Energy Action Plan (SEAP) as of 2017 [31], therefore, both received zero for the first indicator (Table 3) though Veszprém is currently in the process of creating its Sustainable Energy and Climate Action Plan (SECAP). In the case of Veszprém, a cogeneration plant was recently installed [29] while there are only plans for such a system in Zalaegerszeg (though one district heating system operates using geothermal energy) [30]. Energy savings have been accomplished and are continuously implemented in both cities by renovating public buildings and installing photovoltaic systems, e.g. on the flat roofs of a grammar and primary school in Veszprém and on the Mayor’s office in Zalaegerszeg (the performance of which can be accessed online from the webpage of the city). Nonetheless, there is no building with net zero CO2 emissions, that is why both cities received 1 point for the ‘energy savings in end- usage’ indicator. The difference in size of the cities does not necessitate different types of public transport; both Zalaegerszeg and Veszprém have local bus routes that are operated by the same regional bus company. Further points could have been allocated for tram and subway lines (2 for existing, 1 for planned) and an extra point would have been given to the city with the longest tram/subway network [20]. LED technology is considered an efficient public lighting solution (1 point) and an additional point can be gained if solar energy is used to power armatures. Recent investments were made in both cities to improve the energy efficiency of public lighting after the introduction of the cited strategies. 2.1.3. Renewable Energy Potential and Utilization (D3): The renewable energy potential is highly dependent on the location, topology and geology of the area but the local government can have a strong influence on the utilisation of these resources. While regional data could be gathered for the potentials, national data [32] had to be used for the share of renewable sources in terms of electricity production and biofuel use in transport indicators because there was no reliable local information for Veszprém (Table 4). For Zalaegerszeg biogas from the regional municipal wastewater treatment plant is converted to provide the local buses with liquid fuel. Based on a presentation [33] the tanked volume is known for 2015, therefore, the biofuel utilization in terms of transportation was modified accordingly. To attain an accurate comparison the national value for 2015 was considered in the case of Veszprém. Table 4. The data of the Renewable Energy Potential and Utilization (D3) Indicator Veszprém Zalaegerszeg Solar energy potential [Wh/m2/day] 3,425 [25] 3,014 [26] Wind energy potential [m/s] 4.921 [34] 3.505 [34] Geothermal energy potential [mW/m2] 60 [35] 90 [35] Renewable energy usage for electricity [%] 8.76 [32] 8.76 [32] Biofuel utilization in terms of transport [%] 4.15 [32] 5.46 [32-33] Table 2. The data of Energy Consumption and Climate (D1) Indicator Veszprém Zalaegerszeg Energy consumption of buildings [MWh] 309,393 [25] 313,434 [26] Energy consumption of transport (MWh] 225,479 [25] 217,228 [26] Total energy consumption [MWh/capita] (in brackets: population) 19.50 (61,721) [25] 11.65 (59,499) [26] Number of Heating Degree Days (HDD) 2,890 [28] 2,850 [28] Number of Cooling Degree Days (CDD) 1,619 [28] 1,607 [28] Table 3. The data of Penetration of Energy and CO2- Saving Measures (D2) Indicator Veszprém Zalaegerszeg Sustainable Energy Action Plan (SEAP) 0 [31] 0 [31] Combined heat and power-based district heating/cooling system 2 [25] 1 [26] Energy savings in end- usage (buildings) 1 [25] 1 [26] Density of public transport network 1 [29] 1 [30] Efficient public lighting armatures 2 1 ADAPTING THE SDEWES INDEX 45(1) pp. 49–59 (2017) DOI: 10.1515/hjic-2017-0049 53 2.1.4. Water and Environmental Quality (D4): There were no available local data for the domestic blue water footprint, ecological footprint and biocapacity, therefore, national values were applied in the calculation. Air quality is only described in terms of the PM10 concentration [20]. The water quality index was ambiguous as the articles [21-22] referred to the indicator as drinking water quality but the Water Quality Index (WATQI) refers to natural water quality [36]. The index relies on the global database of the United Nations GEMS/Water Programme and includes five indicative parameters: dissolved oxygen, pH, conductivity, total nitrogen and total phosphorous. Unfortunately the WATQI of countries are only available for 2008 [37], from 2012 the water quality index was replaced with access to sanitation and drinking water in terms of the aggregated Environmental Performance. It has to be noted that in the case of Hungary the drinking water is supplied from underground reservoirs that are only linked indirectly to surface waterbodies while in other countries these serve as direct sources of drinking water. Thus using water quality indices for inland water bodies may be good indicators of safe access to water. 2.1.5. Emissions and Industrial Profile (D5): As in the case of the first dimension the emission values of the commercial buildings and municipal vehicle fleet were unavailable, therefore, these could not be included in the calculation. Information on CO2-intense industries was gathered by going through an online company database [41]. The indicator carbon accreditation of airports became zero for both cities for different reasons: Veszprém has no airport and the one near Zalaegerszeg has no accreditation. Since the first case means no emissions while in the second case the existing emissions are not measured, the purpose of the indicator is not fully achieved. The original aim was to include the emissions of the airports in some way in the SDEWES Index as the SEAPs do not take them into consideration [21]. Table 6. The data of CO2 Emissions and Industrial Profile (D5) Indicator Veszprém Zalaegerszeg CO2 emissions of buildings [t CO2] 87,882 [25] 82,579 [26] CO2 emissions of transport [t CO2] 71,711 [25] 70,282 [26] Average CO2 emissions [t CO2/MWh] 0.393 0.292 Number of CO2 intense industries 4 [41] 4 [41] Carbon Accreditation of Airport [levels] 0 0 [42] 2.1.6. City Planning and Social Welfare (D6): Two indicators need further explanation (Table 7). The prices of public transport were introduced in Ref. [21] instead of the share of public transport in terms of total passenger kilometres [22], the latter not being accessible in all cases. The more a single ticket costs, the less likely people will choose public transport. On the other hand, easy access to public transport should result in positive externalities such as cleaner air and less traffic jams, by and large a more liveable city. The urban form and protected sites indicator is an aggregation of several factors (Table 8): compact city form (whether it is mono- or polycentric), urban green areas and surrounding green corridors are evaluated. To determine the compactness of the cities, the energy consumption of transport compared against population density was chosen, as a compact city can be described as of high population density [47] and because of the short distances cars are less likely to be used. Thus the smallest value received 3 points while the highest received 1. Urban green spaces were examined in Table 5. The data of Water and Environmental Quality (D4). Indicator Veszprém Zalaegerszeg Domestic water consumption [m3/capita] 7 [38] 7 [38] Water quality index [/100] 92 [37] 92 [37] Average air quality PM10 [µg/m3] 23.59 [39] 29.60 [39] Ecological footprint [gha] 2.9 [40] 2.9 [40] Biocapacity [gha] 2 [40] 2 [40] Table 7. The data of City Planning and Social Welfare (D6) Indicator Veszprém Zalaegerszeg Price of public transport ticket [EUR] (1 EUR = 310 HUF) 1.07 [43] 1.10 [43] Urban form and protected sites 1 2 GDP per capita [PPP$ national] 25,068.9 [44] 25,068.9 [44] Inequality adjusted well-being (HPI) 4.3 [45] 4.3 [45] Tertiary education rate (national) [%] 21 [46] 21 [46] SEBESTYÉN, SOMOGYI, SZŐKE, AND UTASI Hungarian Journal of Industry and Chemistry 54 comparison with the Hungarian county seats [48]: 0-30 m2/capita: 1 point, 30-50 m2/capita: 2 points and over 50 m2/capita: 3 points. Green corridors were also assessed on a county basis [49] instead of using the suggested GIS- based method [20]. The categories were determined from 1-3 by only taking the green corridor areas of Hungarian counties into consideration. 2.1.7. R&D, Innovation and Sustainability Policy (D7): Results for the seventh dimension are listed in Table 9. The number of public and private universities yielded an unexpected result for the two cities in question. Universities seated in the town and those where only a faculty is based there were equally counted. If only those universities that are seated in the said city were taken into consideration, Veszprém would have 2 versus 0 in Zalaegerszeg. Additional points were given if the university was listed in the Scimago Institutions Rankings [53]. The energy strategy of Veszprém envisions a 25% CO2 emissions reduction by 2026, the basis being 2007 [25], while Zalaegerszeg aims to achieve a 36% reduction by 2030, compared to 2012 [26]. To facilitate a comparison, goals for 2020 were calculated by linear interpolation. 3. Results After processing the necessary calculations, the SDEWES Indices of both Veszprém and Zalaegerszeg were 1.54. As is clear in Fig.1, the values are integers and, except for one case (D5), Veszprém achieved better or equal results. It also has to be noted that on several occasions the difference between the data was very small. Still, the better city was awarded with 1 and the worse value with zero in the normalisation process. In order to eliminate this problem, a third city was included in the benchmark. Ohrid was chosen as the size of this historical Macedonian town is similar to the other two and all data were available from Kilkis [21]. Also, this was the only city in this comparison that had no SEAP. As an alternative solution the ‘average South East European (SEE) city’ from the same article was included to put the two Hungarian cities to the test to see where they would be in the ranking of the SEE cities of that sample. The inclusion of these two examples changed the order of the cities (Table 10). While the results of both Veszprém and Zalaegerszeg improved, Zalaegerszeg gained more from the inclusion of another city from a different country in the benchmark. The reason for the improvement of the indices is that in several cases the national data had to be included in the calculation and hence, the indicator became 0. While both Veszprém and Zalaegerszeg gained points from increasing the sample size, the accumulated increase was larger for Zalaegerszeg (10.03 compared to 8.12). The arrows show in which direction the indicators changed. In the case of Veszprém, data for D3 and D6 decreased but not significantly, for Zalaegerszeg there was no change in D3 and only a slight decline in D5. The SDEWES Index results are similar, the difference between the highest and lowest values is 0.37. Nonetheless, heterogeneity exists with regard to the individual indicators. Results for each dimension are visualised in Fig.2. Both Hungarian towns performed well concerning D1 (energy consumption and climate), D4 (water and environmental quality) and D7 (R&D, innovation and sustainability policy), while there is room for improvement in the fields of CO2-saving measures and city planning. The average SEE city, on the other hand, possesses lower values regarding energy consumption and environmental and water quality. Table 8. Data for Grading Urban Form and Municipal Management Veszprém Zalaegerszeg Urban form and protected sites 1 2 Compact city form 1 3 monocentric x x polycentric population density [capita/km2] 486.38 580.99 Urban green spaces 1 2 urban park intensity [m2/capita] 21 [29] 34.6 [30] Green corridors 1 1 protected sites x x national park/Ramsar x x Table 9. The data of the R&D, Innovation and Sustainability Policy (D7) Indicator Veszprém Zalaegerszeg R&D and innovation policy orientation 3 [50] 2 [50] National patents in clean technologies 2.5 [51] 2.5 [51] Number of public/private universities (city) 3 [52-53] 5 [52-53] h-index of scientific publications 301 [54] 301 [54] Reduction Target for CO2 Emissions (2020) [%] 18 [25] 16 [26] ADAPTING THE SDEWES INDEX 45(1) pp. 49–59 (2017) 55 Figure 1. Results of the first comparison on radar charts. 4. Conclusion The process of gathering data revealed that both Veszprém and Zalaegerszeg need to collect and measure data related to energy efficiency and other indicators of sustainable development more precisely. Creating a database of detailed information on energy use, CO2 emissions and use of renewable sources which is regularly updated would help to achieve the ambitious goal of becoming a sustainable city within a relatively short timespan. Also, the development of a SEAP or SECAP and becoming a member of the Covenant of Mayors would be advantageous and for which Veszprém has started taking steps. Based on the dimensions of the SDEWES Index, Veszprém needs to improve in terms of D3 and D6. The individual indicators highlight that the energy consumption of public transport could be reduced, based on the example of Zalaegerszeg, and also utilization of renewable energy should be improved. In terms of city planning and social welfare the number and area of urban parks can be increased more easily than that of protected sites. Establishing green areas is included within the urban development strategy of the city [29], so improvements may be expected in terms of this indicator. Also, progress in developing a compact city form is anticipated based on the plans to reform the public transport system and relocate the central bus station to next to the railway station [29]. In the case of Zalaegerszeg, dimensions D2, D3 and D6 are lower. Energy-saving measures could be improved by constructing a cogeneration plant to improve the penetration of district heating, using solar panels in public lighting, and also using the wastewater heat to facilitate the full utilization of biogas as a liquid biofuel [55]. In the case of D3, the potential of solar and wind energy cannot be increased and no information concerning the local use of renewable energy resources in terms of electricity was found. As for D6 the same suggestions as in the case of Veszprém can be made to increase the values of the individual indicators. It has to Table 10. Results of calculating the SDEWES Index Veszprém Zalaegerszeg Ohrid Average SEE city D1 3.53 ↑ 3.93 ↑ 3.78 1.43 D2 2.45 ↑ 0.95 ↑ 2.00 3.50 D3 1.99 ↓ 2.00 - 1.60 3.12 D4 4.17 ↑ 3.43 ↑ 2.45 1.16 D5 2.28 ↑ 2.79↓ 3.00 2.45 D6 0.98 ↓ 1.95 ↑ 4.00 3.65 D7 3.72 ↑ 3.97 ↑ 1.08 1.75 SDEWES 2.74 ↑ 2.83 ↑ 2.71 2.47 SEBESTYÉN, SOMOGYI, SZŐKE, AND UTASI Hungarian Journal of Industry and Chemistry 56 Figure 2. Results of the second comparison on radar charts. be pointed out though that the use of the ticket price for public transport resulted in an unjust outcome: since in Ohrid there is no means of local public transport, this indicator became zero, which could also mean that transportation is free. Therefore, Ohrid will always receive the highest value in terms of the process of normalisation as long as there are no local buses in the city. In terms of the process of evaluating the two Hungarian county seats, the benchmarking method was assessed as well. Without a doubt, the SDEWES Index has its benefits. It uses environmental, economic and social indicators, gives credit to CO2 reduction goals and also considers the possible use of renewable resources. Also, human resources are included presuming that higher education and research and development seek to achieve sustainability. On the other hand, the authors identified some drawbacks, too. The first and the fifth dimensions both focus on energy consumption and CO2 emissions. Since these data strongly correlate with each other, the inclusion of both measures leads to redundancy. Also, these two dimensions are weighed more than the others, therefore, energy-related information outweighs other aspects of sustainability. Furthermore, some of the parameters favour smaller cities over larger ones and vice versa. For example if the absolute values of energy consumption of a small city and a capital are compared, the small city will undoubtedly achieve a better result. Figure 3. The evaluation intervals of SDEWES Index indicators ADAPTING THE SDEWES INDEX 45(1) pp. 49–59 (2017) 57 An example of the opposite would be a city with an accredited airport (ACA 3) as opposed to a town with no airport (small towns do not always have airports). Similarly in a capital, where subway and tram lines are at one’s disposal, the density of the public transport network would be high while it would be uneconomical to have trams in a smaller town where the bus lines are sufficient. To overcome the problem of favouring results of cities of different sizes, using data which is proportional to area or population is suggested wherever possible. Another problem was that the scoring of the qualitative indicators such as the urban form and protected sites was not always clear. If the SDEWES Index is to be used widely then these calculations have to be made transparently and be well documented. Due to the nature of the Min-Max method, small differences may be magnified and large differences may diminish. Also, as the index requires certain data that can only be obtained on a national level, the comparison between cities in the same country is somewhat limited. A solution to this problem may be to include towns from different countries and to choose a range (by including more than two cities) in such a way that provides balanced scales in terms of the indicators. While normalisation facilitates the inclusion of values on different scales, the Min-Max method makes it difficult to compare the results of two sets of cities. At present, the extremes are defined by the individual parameters of the cities chosen to be included in the benchmark (Fig.3). In terms of another comparison with a different batch of cities (that can have a common set as with the previous version) the two evaluation intervals may not be equal (��1 ≠ ��2 since �� ∈ ��2 and �� ∉ ��1). As the two extremities are different in the two benchmarks, the results cannot be compared to each other. In the case of extremities a change in the order might appear as was the case in this paper. Including the average of a different batch of cities (given that the average of any parameter is not equal to that of either city) may resolve the limitations of the Min-Max method but only momentarily. Since the ranking is dynamic and changes as the cities develop, the average SDEWES Index of a previous time period will not provide relevant information with regard to the current situation concerning the sample from which the average city was created and neither on the sample of two cities one wanted to expand. Besides the obvious solution of having a sample size of at least three cities, the authors suggest the following: two artificial sets of parameters should be created to serve as absolute extremes of the SDEWES Index. The worst case scenario is referred to as the ‘horror city’ and the best case scenario is named ‘SDEWES city’ after the Index itself. The legend for Fig.3 is as follows:  Ai, Bi, Ci – The measured/real indicator values of the analysed cities in the first calculation.  Si – The measured/real indicator value of the city to be included in the second calculation.  Ii1 – The evaluation interval, when the new city’s value is between the other indicator values.  Ii2 – The evaluation interval, when the new city’s value falls outside of the other indicator values.  zi – The theoretical minimum value of indicator i, the ‘horror city’.  wi – The theoretical maximum value of indicator i, the ‘SDEWES city’. To resolve this problem with regard to the evaluation intervals changing from time to time, the minimum and maximum values of each indicator must be determined in a way that the examined cities could be included in the evaluation intervals: ��1 = ��2 = ��� (5) �� ∈ ��2 , �� ∈ ��1 ... �� ∈ ��� (6) �� ≤ �� ≤ �� (7) Defining these utopian and negative examples requires careful examination of the indicators. Some parameters are dependent on the geographical location while others need to follow a realistic optimal and unfavourable alternative, for example, the tertiary education rate may be zero in the worst case scenario but it is arguable whether 100% would be favourable from the viewpoint of urban management. Further studies are needed to define the ‘horror’ and ‘SDEWES’ cities of the SDEWES benchmarking method. REFERENCES [1] United Nations Department of Economic and So- cial Affairs Population Division: World Urbaniza- tion Prospects: The 2014 Revision, Highlights. (United Nations, New York), 2014 ISBN 978-92-1- 151517-6 [2] United Nations Department of Economic and So- cial Affairs, Population Division: World Urbaniza- tion Prospects: The 2014 Revision. (United Na- tions, New York), 2015 [3] Global Urban Observatory: Global Urban Indica- tors Database. United Nations Human Settlements Programme (UN - Habitat), 2002 [4] Yamato, N.; Sasaki, K.; Hamada, Y.; Ito, K.; Wong, Y.Y.: Global Power City Index 2015, Sum- mary. (Institute for Urban Strategies - The Mori Memorial Foundation, Tokyo, Japan), 2015 [5] Economist Intelligence Unit: European Green City Index. (Siemens AG, Munich, Germany), 2009 [6] Sovacool, B.K.; Brown, M.A.: Twelve metropoli- tan carbon footprints: A preliminary comparative global assessment, Energy Policy, 2010 38, 4856- 4869 DOI: 10.1016/j.proeng.2017.07.146 SEBESTYÉN, SOMOGYI, SZŐKE, AND UTASI Hungarian Journal of Industry and Chemistry 58 [7] Christopher, J.M.; Kammen, D.M.: A Consump- tion-Based Greenhouse Gas Inventory of San Fran- cisco Bay Area Neighborhoods, Cities and Coun- ties: Prioritizing Climate Action for Different Lo- cations, Bay Area Air Quality Management Dis- trict, UC Berkeley, 2015 escholarship.org/uc/item/2sn7m83z, Accessed: 15th August 2016. [8] López-Ruiz, V.R.; Alfaro-Navarro, J.L.; Nevado- Pena, D.: Knowledge-city index construction: An intellectual capital perspective, Expert Syst. Appl., 2014 41, 5560-5572 DOI: 10.1016/j.eswa.2014.02.007 [9] Toro, J.; Requena, I.; Duarte, O.; Zamorano, M.: A qualitative method proposal to improve environ- mental impact assessment, Environ. Impact Assess. Rev., 2013 43, 9-20 DOI: 10.1016/j.eiar.2013.04.004 [10] Eccleston, C.H.: The EIS Book: Managing and Preparing Environmental Impact Statements. (CRC Press), 2013 ISBN 9781466583641 [11] Canter, L.W.: Environmental impact assessment, second edition. (McGraw-Hill, Inc., USA), 1996 ISBN 978-0071141031 [12] Öberg, C.; Huge-Brodin, M.; Björklund, M.: Ap- plying a network level in environmental impact as- sessments, J. Bus. Res., 2012 65, 247-255 DOI: 10.1016/j.jbusres.2011.05.026 [13] Speizer, F.; Magyar, I.; Enisz, K.: Municipal envi- ronmental-monitoring system, Hung. J. Ind. Chem., 2010 38(1), 63-66 [14] Utasi, A.; Yuzhakova, T.; Sebestyén, V.; Németh, J.; Robu, B.; Rédey, Á.; Lakó, J.; Fráter, T.; Rádu- ly, I.; Ráduly, L.; Popita, G.: Advanced Quantita- tive Environmental Impact Assessment Method, Environ. Eng. Manag. J., 2013 12, 305-310 [15] Shaid, A.: Methodological Limitations of Deter- mining Global Pollution Index as a Tool for Envi- ronmental Impact Assessment and a Proposed Ex- tension, Environment and Ecology Research, 2014 2, 240-247 DOI: 10.13189/eer.2014.020603 [16] Saaty, R.W.: The analytic hierarchy process - what it is and how it is used, Math. Modelling, 1987 9(3), 161-176 DOI: 10.1016/0270-0255(87)90473-8 [17] Ching-Lai, H.; Kwangsun, Y.: Multiple Attribute Decision Making - Methods and Applications: A State-of-the-Art Survey. (Springer-Verlag, Berlin, Heidelberg), 1981 ISBN 978-3-642-48318-9 [18] MacCrimon, K.R.: Decision Marking Among Mul- tiple-Attribute Alternatives: a Survey and Consoli- dated Approach RAND memorandum. (The Rand Corporation, Santa Monica, California), 1968 [19] Herva, M.; Roca, E.: Review of combined ap- proaches and multi-criteria analysis for corporate environmental evaluation, J. Clean. Prod., 2013 39, 355-371 DOI: 10.1016/j.jclepro.2012.07.058 [20] SDEWES Index [Online] www.sdewes.org/sdewes_index.php, Accessed: 23rd October 2017. [21] Kılkış, Ş.: Sustainable development of energy, wa- ter and environment systems index for Southeast European cities, J. Clean. Prod., 2016 130, 222- 234 DOI: 10.1016/j.jclepro.2015.07.121 [22] Kılkış, Ş.: Composite index for benchmarking local energy systems of Mediterranean port cities, Ener- gy, 2015 92, 622-638 DOI: 10.1016/j.energy.2015.06.093 [23] Kılkış, Ş.: Benchmarking South East European Cit- ies with the Sustainable Development of Energy, Water and Environment Systems Index, J. Sustain. Dev. Energy Water Environ. Syst., In Press, 2017 DOI: 10.13044/j.sdewes.d5.0179 [24] Holler, A.; Mader-Hirt, D.I.L. Eds.: Climate Star 2012, The European Award for Local Climate Pro- tection Initiatives. (Land Niederösterreich, Gruppe Raumordnung, Umwelt und Verkehr, St. Pölten, Austria), 2012 old.klimabuendnis.org/fileadmin/inhalte/dokumente/20 12/Readme_ClimateStar2012_en.pdf, Accessed: 14th February 2016. [25] Energy Strategy of Veszprém for 2010-2025 (Veszprém Megyei Jogú Város energetikai straté- giája 2010-2025, in Hungarian), Veszprém, 2011 [26] Zalaegerszeg – Ecocity, Renewable Energy Strate- gy (Zalaegerszeg – Ökováros, Megújuló Energia Stratégia, in Hungarian), Zalaegerszeg, 2014 [27] OECD-JRC, Handbook on Constructing Compo- site Indicators: Methodology and User Guide. (OECD, Paris), 2008 ISBN 978-92-64-04345-9 [28] Degree Days, [Online]. www.degreedays.net/, Ac- cessed: 3rd April 2016. [29] Integrated urban development strategy of Veszprém (Veszprém Megyei Jogú Város integrált településfejlesztési stratégia, in Hungarian), Veszprém, 2014 [30] Integrated urban development strategy of Zalae- gerszeg (Zalaegerszeg Megyei Jogú Város integrált településfejlesztési stratégia 2014-2020, in Hungar- ian), Zalaegerszeg, 2014 [31] Covenant of Mayors for Climate & Energy, [Online]. ww.covenantofmayors.eu/actions/sustainable- energy-actionplans_en.html, Accessed: 20th No- vember 2017 [32] International Energy Agency, [Online]. www.iea.org/statistics/statisticssearch/report/?coun try=HUNGARY=&product=balances, Accessed: 9th November 2017 [33] Böcskei, Zs.: Biomethane production in Zalaeger- szeg (presentation), Workshop on procedures ena- bling better use of the energy content of sewage and sewage sludge, Budapest, 22nd June 2017 [34] National Renewable Energy Laboratory, Solar and Wind Energy Resource Assessment, [Online]. maps.nrel.gov/swera/. Accessed: 26th March 2016 [35] Horváth, F.; Bada, G.; Windhoffer, G.; Csontos, L.; Dombrádi, E.; Dövényi, P.; Fodor, L.; Grener- czy, Gy.; Síkhegyi, F.; Szafián, P.; Székely, B.; Timár, G.; Tóth, L.; Tóth, T.: Atlas of the present- day geodynamics of the Pannonian basin: Eurocon- form maps with explanatory text (A Pannon- medence jelenkori geodinamikájának atlasza: Eu- ro-konform térképsorozat és magyarázó, in Hun- garian), Magyar Geofizika, 2006 47(4), 133-137 ADAPTING THE SDEWES INDEX 45(1) pp. 49–59 (2017) DOI: 10.1515/hjic-2017-0049 59 [36] Carr, G.M.; Rickwood, C.J.: Water quality: Devel- opment of an Index to assess country performance. UNEP GEMS/Water Programme 1-24, 2007 [37] Esty, D.C.; Levy, M.A.; Kim, C.H.; de Sherbinin, A.; Srebotnjak, T.; Mara, V.: 2008 Environmental Performance Index, New Haven: Yale Center for Environmental Law and Policy, 2008 [38] Mekonnen, M.M.; Hoekstra, A.Y.: National water footprint accounts: The green, blue and grey water footprint of production and consumption, Value of Water Research Report, 2011 50, 1-50 [39] Hungarian Air Quality Network (Országos Lé- gszennyezettségi Mérőhálózat, in Hungarian), Min- istry of Agriculture, [Online]. www.levegominoseg.hu, Accessed: 3rd April 2016 [40] Global Footprint Network, [Online]. www.footprintnetwork.org/en/index.php/GFN/pag e/trends/hungary/, Accessed: 28th April 2016 [41] Company database (Cégbank, a cégadatbázis, in Hungarian), [Online]. www.cegbank.hu, Accessed: 11th April 2016 [42] Airport carbon accreditation, [Online]. www.airportcarbonaccredited.org/airport/participa nts/europe.html, Accessed: 6th May 2016 [43] North-west Hungary Transportation Center Ltd. (Északnyugat-magyarország Közlekedési Központ Zrt., in Hungarian), [Online]. www.enykk.hu, Ac- cessed: 3rd April 2016 [44] The World Bank, [Online] da- ta.worldbank.org/indicator/NY.GDP.PCAP.PP.CD, Accessed: 4th April 2016 [45] Abdallah, S.; Michaelson, J.; Shah, S.; Stoll, L.; Marks, N.: The Happy Planet Index 2012 Report - A global index of sustainable well-being, London, 2012 ISBN 978 1 908506 17 7 [46] OECD, Country Note - Education at a Glance 2013: Hungary, 2013 DOI: 10.1787/eag-2017-50-en [47] Benton, E.; Jenks, M.; Williams, K.: The Compact City: A Sustainable Urban Form? (Routhledge), 2003 ISBN: 978-0419213000 [48] Hungarian Central Statistical Office, Statistical Yearbook of Hungary 2013 (Központi Statisztikai Hivatal, Magyar Statisztikai Évkönyv 2013, in Hungarian), Budapest, 2014 [49] Hungarian Central Statistical Office, [Online]. Available: www.ksh.hu/docs/hun/xstadat/xstadat_eves/i_ur00 5.html., Accessed: 28th April 2016 [50] Hungarian Central Statistical Office, Research and development 2014: Statistical Reflections (Kutatás- fejlesztés 2014 Statisztikai tükör, in Hungarian), Budapest, 2015 [51] European Patent Office, [Online]. www.epo.org/, Accessed: 3rd April 2016 [52] Site of National Higher Education Information Centre, [Online]. www.felvi.hu, Accessed: 5th May 2016 [53] Scimago Institutions Ranking, [Online]. http://www.scimagoir.com/webvisibility.php?ranki ngtype=webvisibility&indicator=Website%20Size §or=&country=HUN&display=table&page=2 &year=2008., Accessed: 5th May 2016 [54] Scimago Journal & Country Rank, [Online]. http://www.scimagojr.com/countryrank.php, Ac- cessed: 4th April 2016 [55] Baranyák, Z.; Kontra, J.; Kovács, A.; Havas, M.; Jani, I.; Mayer, Z.; Pálfi, Sz.; Garbai, L.; Kovács, K.: Zalaegerszeg Smart City 2050, Goodwill Con- sulting, 2016 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 61–66 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0009 SURFACE ENERGY HETEROGENEITY PROFILES OF CARBON NANOTUBES WITH A COPOLYMER-MODIFIED SURFACE USING SURFACE ENERGY MAPPING BY INVERSE GAS CHROMATOGRAPHY FRUZSINA GERENCSÉR, 1 NORBERT RIEDER, 1 CSILLA VARGA, 2 JENŐ HANCSÓK, 2 AND ANDRÁS DALLOS 1* 1 Department of Physical Chemistry, University of Pannonia, 10 Egyetem str., Veszprém, H-8200, HUNGARY 2 MOL Department of Hydrocarbon and Coal Processing, University of Pannonia, 10 Egyetem str., Veszprém, H-8200, HUNGARY The effectiveness and quantitative control of the surface transition of multi-walled carbon nanotubes (MWCNTs) was characterized by inverse gas chromatography (iGC). The surface energy profile of carbon nanotubes com- patibilized with an olefin-maleic-anhydride-ester-amide (OMAEA)-type coupling agent was determined by a sur- face energy analyzer (SEA). The surface energetic heterogeneity with energy distributions of dispersive and specific (acid-base) components of the surface energy of the MWCNTs were determined at various surface cov- erages. The results of the surface energy mapping showed that surface treatment significantly reduced the dis- persive surface energy of MWCNTs and increased the specific surface energy. Furthermore, the surface modifi- cation enhanced its Lewis basic character and simultaneously decreased the acidic character of MWCNTs. It has been demonstrated that the surface treatment modified the heterogeneity profiles of the energetic surface of the carbonaceous nanomaterials. Keywords: carbon nanotubes, surface treatment, inverse gas chromatography, surface energy analysis 1. Introduction Carbon nanotubes (CNTs) can serve as excellent candi- date materials for uses in numerous industrial applica- tions because of their considerable advantages. CNTs are one of the best reinforcing constituents for nano- composites [1] and hopefully catalytic metal–support in heterogeneous catalysis [2]. CNTs could replace the common catalyst supports of Ni/Mo-catalysts used in the production of fuel components of engine fuels with high hydrogen contents in their molecular structures [3]. Carbon nanotube-supported Co/Mo-catalysts with different Co/Mo atomic ratios were successfully used in the hydrocracking reaction of the vacuum resi- due of crude oil from Gudao oil field [4]. In the Fisch- er–Tropsch process (FTP), CNTs that supported transi- tion metal catalysts are used to increase catalytic activi- ties. An excellent study of FTP on Co catalysts support- ed by CNTs was reported by Tavasoli et al. [5]. Chen et al. [6] demonstrated that Fe nanoparticles encapsulated in CNTs are promising catalyst in FTP to synthesize light olefins. The catalytic consequence of hydrothermal liquefaction of microalgae to produce bio-oil over CNT- *Correspondence: dallos@almos.vein.hu supported transition metal (Co, Ni, Pt) catalysts was reported by Chen et al. [7]. To change the wettability and chemical character of the CNTs or to avoid agglomeration in nanocompo- sites, the CNT surfaces are often exposed to surface functionalization [2] and modification processes using polyfunctional anchoring, capping, and coupling agents [8]. Research has shown that metal–support bindings can be strengthened by functional groups that are cova- lently bonded (grafted) to the support. Functionalized carbon nanotube-supported Pt nanoparticles were ap- plied with favourable results in terms of selective olefin hydrogenation [9]. Because CNTs adsorb molecules well, functionalized CNTs are attractive chromatograph- ic stationary phases for separation of normal and isoal- kanes and aromatic compounds in the development of alternative fuels with high hydrogen/carbon ratios [10]. However, non-covalent functionalization using coupling agents or compatibilizers does not perturb the structure of the carbon nanotubes, establishes proper interactions between carbon nanotubes and the polymer matrix, and prevents the formation of nanotube agglom- erates [11]. In terms of the properties of the reinforced composites of CNTs, the couplings between the nano- tubes and the matrix are important beside the mechani- cal properties of the building parts [12]. These interac- tions depend on the surface properties and energies of the two materials. The surfaces of chemically derivati- zated CNTs were investigated by means of various ana- GERENCSÉR, RIEDER, VARGA, HANCSÓK AND DALLOS Hungarian Journal of Industry and Chemistry 62 lytical methods, e.g. thermal analysis [9-10,13-14], in- frared spectroscopy (IR) [10,13], transmission electron microscopy (TEM) [13,16], Raman [16] and atomic force microscopy [10], and inverse gas chromatography (iGC). iGC is a precise analytical method which is suita- ble for determining the surface energetic characteristics of the CNTs [13,15-17]. iGC was used for the character- ization of the chemical character of the surface and was utilized to measure dispersive and specific surface ener- gies, of numerous CNT substances [18]. The quantita- tive characterization of surface functionalization by surface energy mapping is of great importance. Howev- er, previous papers have presented surface energy val- ues for functionalized CNTs over unclear surface cover- ages without energetic profiles and surface energy dis- tribution functions, which, therefore, could not give correct information on the surface of the CNTs. In this study, the dispersive, specific (acid-base) components of the surface energy with their heterogene- ity charts and energy probability density functions of untreated and compatibilized MWCNTs are presented. A comparative quantitative characterization of the effec- tiveness and quantitative control of surface treatment is given. The exclusive energy scaling of the surfaces of the MWCNTs by energy heterogeneity charts with sur- face energy probability density functions over wide sur- face coverages is the new approach and main novelty of this paper. 2. Experimental and Methods 2.1. Samples and Measurements Multi-walled carbon nanotubes (MWCNTs) were manu- factured at 973 K by the chemical vapour deposition (CVD) process over a Fe/Co bimetallic catalyst at the Department of Chemical Engineering Science (Univer- sity of Pannonia, Veszprém, Hungary) [19]. Their di- ameter was between 10 and 20 nm and their average length was above 30 μm. An olefin-maleic-anhydride-ester-amide (Fig.1) copolymer (OMAEA) was used as a compatibilizer. The coupling agent was synthesized at the Department of MOL Hydrocarbon and Coal Processing (University of Pannonia, Veszprém, Hungary). The surface of MWCNTs was covered by the compatibilizer from a hydrocarbon solution of the coupling agent while the mixture was stirred for 1 hour at 333 K. The solvent was subsequently evaporated and the treated MWCNTs were dried at 383 K for 2 hours in air [11]. The surface energies of as-received and compati- bilized samples of MWCNTs were measured by a Sur- face Energy Analyzer (iGC-SEA, Surface Measurement Systems Ltd., Alperton, UK) over a series of surface coverages from (n/nm) = 0.005 to (n/nm) = 0.030. iGC samples were produced by filling 20-25 mg of CNTs into silanized Pyrex glass tubes of I.D. = 3 mm under a vacuum and moderate vibration. The samples of MWCNTs were stabilized in the column with plugs of silanized glass wool. The samples were preconditioned in the column at the actually measured temperature for 60 minutes before each measurement. The iGC experi- ments were carried out at a column temperature of 353 K, with a Helium carrier gas flow of 10 cm3/min. Me- thane gas was used as a dead-time marker using a flame ionization detector; and n-hexane, n-heptane, n-octane, n-nonane, chloroform and toluene as test compounds. The surface energy values were estimated using the specific retention volumes of the test compounds [20]. The specific retention volumes were obtained from the adjusted retention times: spw mVtV /c ' R  (1) The mean flow rate of the carrier gas in the column, cV , was evaluated as given in Ref. [20]. 2.2. Theoretical Methodologies The dispersive component of the surface energy ( d s ) and its heterogeneity profile of samples of MWCNTs were calculated using the Dorris-Gray method [21] over different surface coverages:   2 CH C,1C, CH d s 22 /ln 4 1             aN VVRT nwnw   (2) When plotting RTln(Vw,nC) against carbon number, nC, for the n-alkane probes, a straight line is generated from the gradient from which the dispersive free energy of the sample surfaces of the MWCNTs, d s , can be calcu- lated. The specific (Lewis acid-base) surface energy ab s of samples of MWCNTs was calculated from the basic component (  s ) and the acidic component (  s ) of the surface energy:  ss ab s 2  (3) The basic and acidic components of the surface energy were obtained from the specific parts of free enthalpy changes of adsorption ab ,ads iG of polar probes i: Figure 1. Structure of the olefin-maleic-anhydride- ester-amide copolymer (OMAEA) coupling agent, where R1: alkyl chain with length of the olefinic mon- omer; R2: alkyl chain with R1–2 carbon number; a, b: 2–21; k: 0.2–2; l: 1–7; m: 1–7 and n: 0.3–2 [11] SURFACE ENERGY HETEROGENEITY PROFILES OF CARBON NANOTUBES 45(1) pp. 61–66 (2017) 63        ss ab ,ads 2  iiii aNG (4) applying the van Oss-Chaudhury-Good theory [22] with the Della Volpe scale [23]. The specific free energy changes of adsorption of the polar probes were obtained as suggested by Donnet et al. [24]. 3. Results and Analysis 3.1. Experiments The dispersive surface energy profiles of the untreated samples of MWCNTs and those treated with the olefin- maleic-anhydride-ester-amide copolymer (OMAEA) coupling agent compatibilized at 353 K and over low surface coverages (n/nm) are presented in Fig.2. The energy profiles show reasonable devaluation of disper- sive surface energy for MWCNTs after surface modifi- cation detected by n-alkane molecular probes: the dis- persive surface energy ( d s ) of the MWCNTs decreased to half of its initial value. The untreated samples of MWCNTs exhibited dis- persive surface energies of ~110 mJ/m2 at 353 K, which is comparable to the values reported by other research- ers studying carbon nanotubes [13-15] and graphitic carbon materials [25]. The relatively high values of the dispersive surface energy of untreated MWCNTs can be attributed to a strong nonpolar interaction potential to build physical long-range Keesom, Debye, and London attractions, which explains their high tendency to ag- glomerate [14]. However, the anchoring of olefin- maleic-anhydride-ester-amide (OMAEA) up on the MWCNTs surface caused a marked decrement in dis- persive part of surface energy from ~110 mJ/m2 to ~48 mJ/m2 at 353 K. The large drop in the value of d s of surface-treated MWCNTs shows that the dispersive surface energy of MWCNTs has been obviously altered by the coupling agent. The surface treatment also affected the dispersive surface energy heterogeneity profile of the MWCNTs. The surface energy mapping of the samples of MWCNTs indicated that the dispersive components of surface energies of untreated samples of MWCNTs are almost constant within the region of low surface cover- age. Consequently, the surface of the untreated MWCNTs can be considered quasi-homogeneous. However, the dispersive surface energy heterogeneity profiles of the treated MWCNTs prove that the copoly- mer-modified MWCNT surface is energetically slightly heterogeneous, because the dependence of d s on sur- face coverage is relatively strong within the region of low surface coverage. In addition, the distributions of the dispersive sur- face energies (Fig.3) obtained by point-by-point integra- tion of dispersive surface energy profiles over the inves- tigated range of the surface coverage support in a more illustrative manner also results in an increase in the dis- persive surface energy heterogeneity. The dispersive surface energy probability function of the MWCNTs became more spread out after modification of the sur- face indicated a greater degree of energetic surface in- homogeneity. The specific surface energy ( ab s ) profiles of the untreated samples of MWCNTs and those treated with OMAEA compatibilized at 353 K and over various sur- face coverages (n/nm) are presented in Fig. 4. The un- treated samples of MWCNTs possess specific surface energy of ~10 mJ/m2 at 353 K, which value is near to that given by Lou et al. (8.84 mJ/m2) for pristine carbon nanotubes at 373 K and over undefined degrees of sur- face coverage [14]. The quantitative surface energy analysis obtained by iGC-SEA methodology demon- strated that surface treatment of MWCNTs resulted in Figure 3. Dispersive surface energy probability func- tions of untreated samples of MWCNTs and those treated with an olefin-maleic-anhydride-ester-amide copolymer (OMAEA) coupling agent compatibilized at 353 K (the solid correlation lines are only to im- prove visualization). Figure 2. Dispersive surface energy profiles of un- treated samples of MWCNTs and those treated with an olefin-maleic-anhydride-ester-amide copolymer (OMAEA) coupling agent compatibilized at 353 K and over various surface coverages (the dotted correla- tion lines are only to improve visualization). GERENCSÉR, RIEDER, VARGA, HANCSÓK AND DALLOS Hungarian Journal of Industry and Chemistry 64 significant changes in surface energies: the specific sur- face energy of CNT surfaces increased more than four- fold, from ~10 mJ/m2 to ~41 mJ/m2. Furthermore, the dependence of ab s on surface coverage is pronounced for the compatibilized samples of MWCNTs which indicates that energetic heteroge- neity attributed to chemical heterogeneity and the exist- ence of electron donor-acceptor atomic groups on the surface. However, the quasi-constant specific compo- nent of surface energy for untreated MWCNTs suggests an energetically homogeneous surface and the absence of high specific energy surface sites. The specific surface energy values of compatibil- ized MWCNTs are much higher than those of the un- treated surfaces, declared the enhanced connection be- tween compatibilized MWCNTs and polar analytes. The observed diversity in specific surface energies is ac- complished from the adsorbed polar atomic clusters: namely the moderately electron-withdrawing maleic anhydride groups; and the electron-donating ester and amide groups with different nucleophilic or electro- philic characteristics. The specific surface energy distributions in Fig. 5 represent the heterogeneity of the samples of MWCNTs and reveal that the untreated MWCNTs exhibited ab s values which varied from 9.9 to 10.3 mJ/m2. This small variation in the specific surface energy demonstrated a fairly energetically homogeneous surface for untreated MWCNTs. However, as surface treatment increased the concentrations of polar clusters on the surface of MWCNTs, the modified surface exhibited great varia- tions in ab s (from 37.3 to 42.9 mJ/m2), implying that the compatibilized MWCNTs are surface energetic het- erogeneous. The surface treatment also modified the chemical characteristics of the MWCNTs. The acid-base surface energy mapping of the samples of MWCNTs (Figs.6 and 7) indicate that both the electron-accepting and do- nating abilities of the MWCNTs were raised apprecia- bly after compatibilization of untreated MWCNTs con- sistent with the adsorption of electron-withdrawing and electron-donating atomic groups on the surface. The untreated samples of MWCNTs possess base surface energy of ~16 mJ/m2 at 353 K, which value is near to that given by Lou et al. (12.97 mJ/m2) for pristine car- bon nanotubes at 373 K and over undefined surface coverage [14]. The surface energy mapping using the iGC-SEA methodology confirms that the surface treatment of MWCNTs raised the basic component (  s ) of surface energy of MWCNTs: it resulted in a sevenfold increase from ~16 mJ/m2 to ~112 mJ/m2. The acidic component of surface energy  s of untreated samples of MWCNTs was measured as ~1.7 mJ/m2 at 353 K, which is in good agreement with that reported by Lou et al. (1.51 mJ/m2) for pristine carbon nanotubes at 373 K and over unde- fined surface coverages [14]. The surface modification of MWCNTs resulted in a more than twofold increase in the value of the acidic component (  s ) of the surface energy of MWCNTs from ~1.7 mJ/m2 to ~3.7 mJ/m2. However, the profiles of acid-base surface energy heterogeneity of the treated MWCNTs also prove that the copolymer-modified MWCNT surface became en- ergetically more heterogeneous, because the dependen- cies of  s and  s on surface coverage are relatively stronger than those of the untreated MWCNTs. The presence of non equi-energetic active surface centers exposes that the surfaces of the compatibilized MWCNTs are not energetically uniform for specific acid-base interactions and the surface treatment consid- erably altered the ability of MWCNTs to connect with molecular species by specific interactions. Figure 5. Specific surface energy probability functions of untreated samples of MWCNTs and those treated with an olefin-maleic-anhydride-ester-amide copoly- mer (OMAEA) coupling agent compatibilized at 353 K (the solid correlation lines only improve visualiza- tion). Figure 4. Specific surface energy profiles of untreated samples of MWCNTs and those treated with an olefin- maleic-anhydride-ester-amide copolymer (OMAEA) coupling agent compatibilized at 353 K and over vari- ous surface coverages (the dotted correlation lines on- ly improve visualization). SURFACE ENERGY HETEROGENEITY PROFILES OF CARBON NANOTUBES 45(1) pp. 61–66 (2017) 65 The larger change in the base component of the surface energy (  s ) after compatibilization with a pol- yalkenyl-poly-maleic-anhydride-ester-amide additive indicates a higher concentration of electron-donating ester and amide groups on the surface. The large value of  s (relative to  s ) implies a more basic characteris- tic and donor properties of the surface of the treated MWCNTs. 4. Conclusion The experimental data demonstrated that iGC is a useful methodology of characterizing the variation in the sur- face characteristics of MWCNTs after non-covalent functionalization. The exclusive energy scaling of the SEA methodology using energy heterogeneity charts with surface energy probability density functions over wide surface coverages presents profitable additional information on the differences in terms of the nature, homogeneity and heterogeneity of surface energies re- sulting from surface transformations. The multilateral surface energy analysis of SEA presents a quantitative control of the effectiveness of surface treatment and demonstrates the importance of the dependence of sur- face energy analysis on coverage. Acknowledgement The present work was published within the framework of the project GINOP-2.3.2-15-2016-00053. SYMBOLS 2CHa , ai cross sectional area of an adsorbed methylene group and of probe i ab ,ads iG specific free enthalpy of adsorption msp mass of the adsorbent in the column N Avogadro’s number R gas constant T temperature ' Rt adjusted retention time cV mean flow rate of the carrier gas Vw specific retention volume ab s , d s specific and dispersive parts of the surface energy of solid sample materi- al 2CH surface energy of a methylene group  s ,  s acid-base components of the surface energy of solid sample material  i ,  i acid-base components of the surface tension of polar liquid probe i Θ=n/nm surface coverage REFERENCES [1] Dresselhaus, M.S.; Dresselhaus, G.; Avouris, P.: Carbon nanotubes: synthesis, structure, properties, and applications (Springer-Verlag GmbH, Berlin, Germany) 2001 DOI: 10.1007/3-540-39947-X Figure 6. Basic component profiles of surface energy of untreated samples of MWCNTs and those treated with an olefin-maleic-anhydride-ester-amide copoly- mer (OMAEA) coupling agent compatibilized at 353 K and over various surface coverages (the dotted cor- relation lines only improve visualization). Figure 7. Acidic component profiles of surface energy of untreated samples of MWCNTs and those treated with an olefin-maleic-anhydride-ester-amide copoly- mer (OMAEA) coupling agent compatibilized at 353 K and over various surface coverages (the dotted cor- relation lines only improve visualization). GERENCSÉR, RIEDER, VARGA, HANCSÓK AND DALLOS Hungarian Journal of Industry and Chemistry 66 [2] Yan, Y.; Miao, J.; Yang, Z.; Xiao, F.X.; Yang, H.B.; Liu, B.; Yang, Y.: Carbon nanotube cata- lysts: recent advances in synthesis, characterization and applications, Chem. Soc. Rev., 2015 44, 3295- 3346 DOI: 10.1039/C4CS00492B [3] Sági, D.; Holló, A.; Varga, G.; Hancsók, J.: Co- hydrogenation of fatty acid by-products and differ- ent gas oil fractions, Journal of Cleaner Produc- tion, 2017 161, 1352-1359 DOI: 10.1016/j.jclepro.2017.05.081 [4] Li, C.; Shi, B.; Cui, M.; Shang, H.J.; Que, G.H.: Application of Co-Mo/CNT catalyst in hydro- cracking of Gudao vacuum residue, J. Fuel. Chem. Technol., 2007 35(4), 407−411 DOI: 10.1016/S1872-5813(07)60026-7 [5] Tavasoli, A.; Sadagiani, K.; Khorashe, F.; Seifkor- di, A.A.; Rohani, A.A.; Nakhaeipour, A.: Cobalt supported on carbon nanotubes — A promising novel Fischer–Tropsch synthesis catalyst, Fuel Process. Technol., 2008 89(5), 491–498 DOI: 10.1016/j.fuproc.2007.09.008 [6] Chen, X.; Deng, D.; Pan, X.; Bao, X.: Iron catalyst encapsulated in carbon nanotubes for CO hydro- genation to light olefins, Chinese J. Catal., 2015 36(9), 1631–1637 DOI: 10.1016/S1872-2067(15)60882-8 [7] Chen, Y.; Mu, R.; Yang, M.; Fang, L.; Wu, Y.; Wu, K.; Liu, Y.; Gong, J.: Catalytic hydrothermal liquefaction for bio-oil production over CNTs sup- ported metal catalysts, Chem. Eng. Sci., 2017 161, 299–307 DOI: 10.1016/j.ces.2016.12.010 [8] Varga, Cs.; Bartha, L.: A novel route for injection moulding of long carbon fibre-reinforced LLDPE, J. Reinf. Plast. Comp., 2014 33(20), 1902-1910 DOI: 10.1177/0731684414549443 [9] Chen, P.; Chew, L.M.; Xia, W.: The influence of the residual growth catalyst in functionalized car- bon nanotubes on supported Pt nanoparticles ap- plied in selective olefin hydrogenation, J. Catal., 2013 307, 84–93 DOI: 10.1016/j.jcat.2013.06.030 [10] Speltini, A.; Merli, D.; Quartarone, E.; Profumo, A.: Separation of alkanes and aromatic compounds by packed column gas chromatography using func- tionalized multi-walled carbon nanotubes as sta- tionary phases, J. Chromatogr. A, 2010 1217(17), 2918–2924 DOI: 10.1016/j.chroma.2010.02.052 [11] Szentes, A.; Varga, Cs.; Horváth, G.; Bartha, L.; Kónya, Z.; Haspel, H.; Szél, J.; Kukovecz, Á.: Electrical resistivity and thermal properties of compatibilized multi-walled carbon nano- tube/polypropylene composites, Express Polym. Lett., 2012 6(6), 494–502 DOI: 10.3144/expresspolymlett.2012.52 [12] Varga, Cs.; Miskolczi, N.; Bartha, L.; Lipóczi, G.; Falussy, L.: Improving the compatibility of man- made fibre reinforced composites, Hun. J. Ind. Chem., 2008 36(1-2), 137-142 [13] Zhang, X.; Yang, D.; Xu, P.; Wang, C.; Du, Q.: Characterizing the surface properties of carbon nanotubes by inverse gas chromatography, J. Ma- ter. Sci., 2007 42(17), 7069–7075 DOI: 10.1007/s10853-007-1536-7 [14] Luo, Y.; Zhao, Y.; Cai, J.; Duan, Y.; Du, S.: Effect of amino-functionalization on the interfacial adhe- sion of multi-walled carbon nanotubes/epoxy nanocomposites, Materials and Design, 2012 33, 405–412 DOI: 10.1016/j.matdes.2011.04.033 [15] Menzel, R.; Tran, M.Q.; Menner, A.; Kay, C.W.M.; Bismarck, A.; Shaffer, M.S.P.: A versa- tile, solvent-free methodology for the functionali- sation of carbon nanotubes, Chem. Sci., 2010 1(5), 603-608 DOI: 10.1039/C0SC00287A [16] Menzel, R.; Lee, A.; Bismarck, A.; Shaffer, M.S.P.: Inverse gas chromatography of as-received and modified carbon nanotubes, Langmuir, 2009 25(14), 8340–8348 DOI: 10.1021/la9000607s [17] Díaz, E.; Ordóñez, S.; Vega, A.: Characterization of nanocarbons (nanotubes and nanofibers) by in- verse gas chromatography, J. Phys. Conf. Ser., 2017 61(1), 904-908 DOI: 10.1088/1742-6596/61/1/180 [18] Menzel, R.; Lee, A.; Bismarck, A.; Shaffer, M.S.P.: Deconvolution of the structural and chem- ical surface properties of carbon nanotubes by in- verse gas chromatography, Carbon, 2012 50(10), 3416–3421 DOI: 10.1016/j.carbon.2012.02.094 [19] Szentes, A.; Horváth, G.: Role of catalyst support in the growth of multi-walled carbon nanotubes, Hun. J. Ind. Chem., 2008 36(1-2), 113–117 [20] Kondor, A.; Dallos, A.: Adsorption isotherms of some alkyl aromatic hydrocarbons and surface en- ergies on partially dealuminated Y faujasite zeolite by inverse gas chromatography, J. Chromatogr. A, 2014 1362, 250-261 DOI: 10.1016/j.chroma.2014.08.047 [21] Dorris, G.M.; Gray, D.G.: Adsorption of n-alkanes at zero surface coverage on cellulose paper and wood fibers, J. Colloid Interf. Sci., 1980 77(2), 353-362 DOI: 10.1016/0021-9797(80)90304-5 [22] van Oss, C.J.; Good, R.; Chaudhury, M.K.: Addi- tive and nonadditive surface tension components and the interpretation of contact angles, Langmuir, 1988 4(4), 884-891 DOI: 10.1021/la00082a018 [23] Della Volpe, C.; Sibioni, S.: Some reflection on acid-base solid surface free energy theories, J. Col- loid Interf. Sci., 1997 195(1), 121-136 DOI: 10.1006/jcis.1997.5124 [24] Donnet, J.B.; Park, S.J.; Balard, H.: Evaluation of specific interactions of solid surfaces by inverse gas chromatography. A new approach based on po- larizability of the probes, Chromatographia, 1991 31, 435-440 DOI: 10.1007/BF02262385 [25] Donnet, J.B.; Park, S.J.: Surface characteristics of pitch-based carbon fibers by inverse gas chroma- tography method, Carbon, 1991 29(7), 955-961 DOI: 10.1016/0008-6223(91)90174-H HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 67–71 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0010 INITIAL ELECTRICAL PARAMETER VALIDATION IN LEAD-ACID BATTERY MODEL USED FOR STATE ESTIMATION BENCE CSOMÓS, DÉNES FODOR, * AND GÁBOR KOHLRUSZ Department of Automotive Mechatronics, Institute of Mechanical Engineering, University of Pannonia, Egyetem u. 10., Veszprém, H-8200, HUNGARY The paper presents a current impulse-based excitation method for lead-acid batteries in order to define the ini- tial electrical parameters for model-based online estimators. The presented technique has the capability to track the SoC (State of Charge) of a battery, however, it is not intended to be used for online SoC estimations. The method is based on the battery’s electrical equivalent Randles’ model [1]. Load current impulse excitation was applied to the battery clamps during discharge while the voltage and current was logged. Based on the Randles’ model, a model function and a fit function were implemented and used by exponential regression based on the measured data. The diffusion-related non-linear characteristic of the battery was approximated by a capacitor- like linear voltage function for speed and simplicity. The initial capacitance of this bulk capacitor was estimated by linear regression on measurements recorded in the laboratory. Then, the RC parameters of the equivalent battery model were derived from exponential regression on transients during each current impulse cycle. The battery model with initial RC parameters is suitable for model-based online observers. Keywords: Battery, SoC, Exponential regression, Randles’ model, Load current impulse 1. Introduction In our daily lives, the number of mobile devices and utilities that can operate without grid connections is increasing. Even though lithium batteries possess better performance properties and energy indicators, lead-acid batteries are still cheaper, significantly present in commercial applications and almost fully recyclable. Therefore, any developments in lead-acid battery systems are still of interest. According to Ref. [1], several methods exist to estimate a battery’s State of Charge (SoC) and State of Health (SoH) but model-based prediction is the most widespread because of its reliablity and robustness. Model-based methods, as the name suggests, need a valid, properly detailed electric battery model. The Randles’ model as a standard battery model is very popular in the contexts of lead-acid and lithium-ion batteries because of its cost-effectiveness and the similarities of both types. By similarity it is meant that the same model can be reasonably used for the parameter estimation of both battery types [2-3]. Some additions to the standard Randles’ model can be made if more details in electrochemistry are required such as diffusion in the bulk and porosity amongst others. The model requires values of initial resistance (R) and capacitance (C). The more accurate the initial *Correspondence: fodor@almos.vein.hu parameters of the model, the faster and more reliable the convergence of a model-based predictor to the actual state, that is, the actual SoC. The scope of the present work is to identify the initial values of RC components (parameters) by evaluating the voltage impulse responses excited by load currents in the time domain. 2. Battery model for impulse excitation In this paper, a standard Randles’ model [7-12] was analyzed that consists of charge-transfer resistance, Rct, battery serial resistance, Rs, double-layer capacitance, Cdl, and bulk capacitance, Cb (Fig.1). The voltage references of the capacitors and currents in Fig.1 were set for discharge. Figure 1. Randles’ battery model applied to a discharging battery pack Rct Cdl Cb Udl Ub ibi(t) Uocv(t) Rs ict idl CSOMÓS, FODOR, AND KOHLRUSZ Hungarian Journal of Industry and Chemistry 68 2.1. State-space model and model function By neglecting the intermediate mathematical steps and rearranging the state-space model, the system can be written in the following form: i C C U U RC U U dt d                                   b dl b dl ctdl b dl 1 1 00 0 1 (1) iR U U U s10 01 ocv b dl              . (2) By solving the output in Eq.(2) in terms of the time domain using a current impulse as the system input, the output function in Eq.(2) can be expressed by Eq.(3) that can be considered as the model function of the system. This form of the output equation serves as a basis for creating a fit function on measured voltage data and then, derive R and C values from the fit parameters. For clarity, each of the terms of the output equation are grouped by alphabetical letters: DBtAetU t   )(ocv (3) where tags A, B and C can be written according to )(ctdl0 tiRUA  (4) t U C ti B 0b b )(  (5) )()( cts tiRRD  (6) where Uocv is the battery’s open-circuit voltage,  is the system’s time constant, t is the measurement time, and i(t) is the impulse load current. Since the proposed method is based on load and relaxation cycles that follow each other during the analysis, Ub0 represents the initial voltage of the bulk capacitor while Udl0 is the initial voltage of the double-layer capacitor at the beginning of each impulse cycle. It should be noted that changes in current during each impulse cycle can be neglected as a result of working in the short time-constant region of the discharge curve, thus the current can be considered as a constant. The battery model shown in Fig.1 is prepared for short-time transient analysis. Even though the model can be used and remains valid for modelling discharge processes that last for several hours, that is, for long- time transients, the accuracy of the model becomes poor under such circumstances. The reason for this is that the battery model presented excludes the diffusion effect that can even be observed by the initial valley-like voltage drop and later as a circle-like voltage response on the long-time discharge voltage curve (Fig.2). Such an exclusion was made because the scope of the current work is to focus on short-time voltage responses to avoid excessive measurement time intervals and computational resources. Consequently, the battery model was optimized for fast SoC detection by short- time battery checks. 2.2. Determining the initial capacity of Cb In Eq.(1), the Ub voltage represents the equilibrium voltage of the battery, therefore, it is related to the battery’s main charge and as a result its SoC. If the battery is excited by small C/10 - C/30 currents, Ub and Uocv can be considered equal. Instead of Ub, Uocv can be measured at the battery terminals. Therefore, the relationship between Uocv and SoC is important and can be determined from laboratory OCV measurements. A small discharge current was applied to the battery terminals under controlled conditions until the battery’s OCV reached the factory’s minimum voltage threshold from a fully charged state. The voltage, current and temperature were logged while the SoC could be calculated by the simple Coulomb-counting method during the process and saved in a lookup table. Then, the lookup table could be used to determine a discrete relationship betwen Uocv and SoC in itself. The Uocv - SoC characterisation method can be extended by the regression on lookup data in order to create a continous Uocv - SoC relationship. A linear function, such as a capacitor-like regression of Uocv - SoC characteristics can be legitimate if the battery’s excitation current is small, i.e. between C/10 and C/15 and is not discharged under 20-25% of SoC. This could be the case when low-power devices are considered as loads. In the case of plain discharge, the SoC changes can be basically tracked by the basis of the B term like in Eq.(5) as conducted in the Coulomb-counting method [1]. In the presented battery model, the B term provides information on the long-term state of the battery and requires initial parameters such as Cb and Ub0. The former represents the battery’s initial capacity, the latter is related to the battery’s initial voltage at the beginning of each impulse cycle. Right before the first current impulse, Ub0 is equal to the battery’s equilibrium voltage since a relaxation time of between 30 minutes and one hour is sufficient for chemical processes to decay. Figure 2. Discharge characteristics of a 15Ah AGM battery excited by different load currents INITIAL ELECTRICAL PARAMETER VALIDATION IN LEAD-ACID BATTERY MODEL 45(1) pp. 67–71 (2017) 69 The initial value of Cb is crucial because it has a strong influence on both the initial voltage drop and the gradient of the long-time discharge voltage curve. In the proposed model, a simplified, that is, capacitor-like formula was used for initial battery capacity estimation so it approximates the battery non-linear discharge characteristic linearly. The introduction of a regression error of a few percent, however, can lead to an easier and faster determination of the initial capacity Cb. The calculation of the battery’s initial capacity was realized according to [4]. The fully charged battery at room temperature needed to be slowly discharged by a C/15 current until its OCV voltage reached the factory recommended minimum voltage threshold. Once the discharge had finished, 2 hours of relaxation had to be observed in order to return the battery back to its almost equilibrium state. Then, a C/15 slow charge had to proceed until the battery’s OCV voltage reached the factory recommended maximum voltage threshold. In both cases, the battery’s OCV voltage, current and temperature were recorded (Fig.3/a). After averaging the discharge and charge-voltage measurements, linear regression was conducted on it according to 0bb b 15/C| )( 1 )( UtQ C tUocv  (7) where Qb(t) can be estimated by Qb(t) = i0 t. Eq.(7) can be identified by the standard form of the linear curve y=mx+b. In Eq.(7), Cb is the battery’s initial capacity, Qb(t) is the actual charge of the battery and Ub0 is the actual voltage of Cb at the beginning of the impulses. By rearranging Eq.(7), Cb can be determined. Even though charge/discharge currents and SoC levels are constrained, the ambient temperature should be controlled as well to provide a constant temperature during the test periods. The value of Cb was calculated as 37766 F at 22°C using C/15 load currents. 3. Exponential regression to derive RC model parameters The evaluation of measurement data and the comparision of measurements and simulations were realized using Matlab. The RC parameters in the model shown in Fig.1 and later in an OrCAD circuit were derived by an exponential regression using an appropriate fit function on the measured voltage data. The regression error between the measured and modelled characteristics can be minimized if the fit function follows the form of the model function, that is, both of them implement similar dynamics and the physical background of the inspected system. Therefore, the fit function can be written as BeAtU t ˆˆ)( ˆ ocv     (8) in a form that is similar to Eq.(3). The hat sign means that the form of Eq.(8) is similar to Eq.(3), but uses a different reference system. When using Eq.(8) attention must be paid to the determination of RC parameters. Eq.(8) gives voltage references with respect to the ground, that is the x-axis, while Eq.(3) yields the references Udl and Us which correspond to Ub. References used by Eqs.(3) and (8) must be matched to derive correct RC parameters (Fig.3/b). It can be seen that Eq.(8) neglects the linear term from Eq.(3). Since the effect of continuous slow discharge of the battery, which is linked to the Cb bulk capacitance in the battery model, possesses a time- constant several orders of magnitude higher than the Cdl-Rct subsystem, it cannot be observed during the short impulse cycles. However, it should be considered during the whole discharge process so Eqs.(3) and (8) need to be used together to estimate the SoC during long-term discharge processes. According to Refs. [5-6], it is practical to evaluate only the discharge component of the voltage responses during every load cycle. This method is also supported by the fact that the dynamic behaviour of the battery for both load and relaxation states can be described by the same battery model since both responses originate from the same battery structure. Because load curves are only needed, the measured voltage should be separated from global data and then, concatenated right after each other. Since voltage data has been prepared in this way, a b Figure 3. (a) Linear regression on the average of C/15 discharge and charge voltage data to derive initial value of Cb. Initial cut-off has been filtered. (b) References of model function and fit function CSOMÓS, FODOR, AND KOHLRUSZ Hungarian Journal of Industry and Chemistry 70 regression can be made during every discharge impulse cycle simultaneously. The relaxation time during each impulse was set to provide enough time for the battery’s double-layer capacitance to almost fully discharge. It is practical because it simplifies Eq.(8) by reducing the effect of the Udl0 term in Eq.(4). The load time was set taking into consideration the time constant of the Cdl - Rct subsystem. According to experiments, it is within the range of 1 - 4s. The calculation method of the term B̂ in (8) is based on the calculation of the limit of the exponential term in Eq.(8). Practical experiments have proven that maintaining a load cycle of 4 in length can speed up and correct the regression. The RC values change during discharge. Due to offline voltage response evaluation, it is not possible to optimise these parameters according to load conditions. Therefore, the average RC values can be used for the whole period of discharge though this introduces a slight misalignment between the measured and modelled voltage characteristics. By applying Eq.(8) on the prepared measurement data, the average RC parameters shown in Table 1 were derived. 4. Model implementation and validation of the impulse excitation method in OrCAD The aim of the OrCAD implementation was to validate the proposed battery model and the applicability of the impulse excitation method for RC parameter determination. The equivalent circuit model shown in Fig.1 was transformed into an electric circuit. The duration of the simulation was set to 1 hour and the time step was equal to the sample time of the real measurement, namely 100 ms. 4.1. Initialization of the OrCAD simulator and RC elements The initial voltages of the capacitors and the proper directions of OrCAD elements should be set carefully. OrCAD uses a reference system of its own that influences the current references of each component. This should be considered while placing an element in the circuit editor and when comparing the electric circuit references to the Randles’ model. The initial voltages of the capacitors, Udl and Ub, were set to describe the battery’s discharge process and thus also follow the voltage references of the Randles’ model, shown in Fig.1. The values of initial capacitor voltages were derived from the chemical background and assumptions. Cdl can be considered fully discharged through Rct after the battery had relaxed (no load) for 2 hours so its initial voltage, Udl, was set to 0. The bulk capacitance of the battery, Cb, reflects the lengthy time constant as well as diffusion-related processes and is related to the main charge storage capability of the battery. If the relaxation time is sufficiently long, the battery reaches an equilibrium state and its OCV becomes equal to Ub. An optimal relaxation time that is sufficiently long was derived from preceding measurements of discharge/charge cycles by applying different load currents and relaxation times. Based on experiments, a relaxation time of 2 hours was applied and as a result the OCV could be considered equal to Ub. Using these assumptions, the initial voltage, Ub, was set to 12.7 V and Udl to 0 V. The introduction of the load current as an impulse excitation can be achieved through a switched power resistor element which is setup in a similar way to the real arrangement (Fig.3/b). The switching routine of the OrCAD element was tuned in accordance with real load-relaxation cycles that were applied to the real testbench. The resistor sets the load current that discharges the battery. During this work, a load current of 3A was used with a load of 5s and relaxation cycles that lasted 10s. All of the initial values are summarized in Table 2. 4.2. Comparing the simulation with the battery model In Fig.4/a, the results of measured and simulated voltage responses are shown. The validity of the model was analyzed by the comparison of measured and simulated voltages. According to the setup, the comparison is performed within a time frame of 5,600 s. The blue curve that represents the simulated data has a longer tail than the red one because filtering needed to be performed on the measured data to cut the initialization process at the beginning of the testbench. The zoomed-in segment shows the fitness of the simulated voltage response. The difference between the curves is relatively small, around 0.05 V to be precise. This error occurred because average RC values were used in the model instead of online tuned ones. Table 2. Initial values of simulation parameters that need to be set before a run Simulation parameter Initial value Simulation step time 100 ms Simulation duration 1 hour Udl0 0 V Ub0 12.7 V Relaxation time / cycles 10 s Load time / cycles 5 s Load current 3 A (with a 4 Ohm load) Table 1. Derived RC values of the battery model from regression Derived parameter name Average value Rct 0.032 Ohm Cdl 92 F Rs 0.056 Ohm SoC after discharge 87.9 % INITIAL ELECTRICAL PARAMETER VALIDATION IN LEAD-ACID BATTERY MODEL 45(1) pp. 67–71 (2017) 71 5. Conclusions This work shows a current impulse-based excitation method that can be used to either track SoC changes during moderate discharge or find proper RC model values for model-based algorithms. The method is founded on an equivalent model-based approach that can implement the dynamic behaviour of a battery without using excessive chemical equations. The technique uses offline analysis of a battery, therefore, it is not able to reasonably track SoC changes in real-time. This technique is not intended to estimate the SoC of batteries by itself, indeed, it estimates good set points for online estimators such as Kalman filters or other model-based observers. The advantages of the presented approach are its rapid nature and simplicity while minimising the error of SoC estimation. The disadvantages of this method are its offline nature and related consequences. This technique could potentially be applied to embedded systems and commercially. 6. Acknowledgement This research was supported by the EFOP-3.6.1-16- 2016-00015 project. The project was supported by the Hungarian Government and co-financed by the European Social Fund. The project was supported by the European Union and co-financed by the European Social Fund (EFOP- 3.6.2-16-2017-00002). REFERENCES [1] Piller, S.; Perrin, M.; Jossen, A.: Methods for state- of-charge determination and their application, Journal of Power Sources, 2001 96(1), 113–120 DOI: 10.1016/S0378-7753(01)00560-2 [2] Spagnol, P.; Rossi, S.; Savaresi, S.M.: Kalman fil- ter SoC estimation for Li-ion batteries, IEEE Inter- national Conference on Control Applications, 28- 30 Sept. 2011, Denver, USA, ISSN 1085-1992 [3] Devarakonda, L.; Wang, H.; Hu, T.: Parameter identification of circuit models for lead-acid batter- ies under non-zero initial conditions, American Control Conference, 4-6 June 2014, Portland, USA [4] Platt, G.L.: Battery modeling, Vol. I., 2015 ISBN-13: 978-1-63081-023-8 [5] Banaei A.; Fahimi B.: Real time condition moni- toring in Li-ion batteries via battery impulse re- sponse, IEEE Vehicle Power and Propulsion Con- ference (VPPC), 1-3 Sept. 2010, Lille, France ISSN 1938-8756 [6] Banaei A.; Khoobroo A.; Fahimi B.: Online detec- tion of terminal voltage in Li-ion batteries via im- pulse response, IEEE Vehicle Power and Propul- sion Conference, 2009. VPPC '09, 7-10 Sept. 2009, Dearborn, USA ISSN 1938-8756 [7] Kularatna, N.: Dynamics and modeling of re- chargeable batteries, IEEE Power Electronics Magazine, 2014 1(4), 23–33 DOI: 10.1109/MPEL.2014.2361264 [8] Lee, Y.-D.; Park, S.-Y.; Han, S.-B.: Online embed- ded impedance measurement using high-power battery charger, IEEE Transactions on Industry Applications, 2015 51(1), 498-508 DOI: 10.1109/TIA.2014.2336979 [9] Kiani, M.: High resolution state of charge estima- tion in electrochemical batteries, Applied Power Electronics Conference and Exposition (APEC), Twenty-Eighth Annual IEEE, 17-21 March 2013, Long Beach, USA [10] Kujundzic, G.; Ileš, Š.; Matuško, J.; Vašak, M.: Optimal charging of valve-regulated lead-acid bat- teries based on model predictive control, Applied Energy, 2017 187 189–202 DOI: 10.1016/j.apenergy.2016.10.139 [11] Bhangu B.S.; Bentley, P.; Stone, D.A.; Bingham C.M.: Observer techniques for estimating the state- of-charge and state-of-health of VRLABs for hy- brid electric vehicles, Vehicle Power and Propul- sion 2005 IEEE Conference, 2005 DOI: 10.1109/VPPC.2005.1554646 [12] Wang, Y.; Fang, H.; Zhou, L.; Wada, T.: Revisit- ing the state-of-charge estimation for lithium-ion batteries, IEEE Control Systems Magazine, 2017 37(4), 73–96 DOI: 10.1109/MCS.2017.2696761 a b Figure 4. (a) Validation of modelled (OrCAD) and measured voltage data in Matlab environment (b) Battery testbench used during the analyzis HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 73-84 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0011 SIMULATING ION TRANSPORT WITH THE NP+LEMC METHOD. APPLICA- TIONS TO ION CHANNELS AND NANOPORES. DÁVID FERTIG ∗1, ESZTER MÁDAI1, MÓNIKA VALISKÓ1, AND DEZSŐ BODA1,2 1Department of Physical Chemistry, University of Pannonia, Egyetem u. 10., Veszprém, H-8200, HUNGARY 2Institute of Advanced Studies Kőszeg (iASK), Chernel u. 14., Kőszeg, H-9730, HUNGARY We describe a hybrid simulation technique that uses the Nernst-Planck (NP) transport equation to compute steady-state ionic flux in a non-equilibrium system and uses the Local Equilibrium Monte Carlo (LEMC) simulation technique to es- tablish the statistical mechanical relation between the two crucial functions present in the NP equation: the concentration and the electrochemical potential profiles (Boda, D., Gillespie, D., J. Chem. Theor. Comput., 2012 8(3), 824–829). The LEMC method is an adaptation of the Grand Canonical Monte Carlo method to a non-equilibrium situation. We apply the resulting NP+LEMC method to ionic systems, where two reservoirs of electrolytes are separated by a membrane that al- lows the diffusion of ions through a nanopore. The nanopore can be natural (as the calcium selective Ryanodine Receptor ion channel) or synthetic (as a rectifying bipolar nanopore). We show results for these two systems and demonstrate the effectiveness of the NP+LEMC technique. Keywords: ion transport, Nernst-Planck, Monte Carlo simulation, nanopore, ion channel 1. Introduction We dedicate this paper to honoring Professor János Liszi, who was the supervisor of one of the authors (DB) and re- spected senior to another (MV), supporting to the careers of both. We (DB and VM) publish this paper together with our students (DF and EM) to demonstrate that the lesson of Professor Liszi — one of the most important goals of senior scientists is to pave the way to the junior ones — has not been left unconsidered. The goal of this work is to present a computer sim- ulation technique for computing steady-state ion trans- port that was developed and applied in our research group with the essential help of Dirk Gillespie [1]. The method is based on the Nernst-Planck (NP) transport equation coupled to the Local Equilibrium Monte Carlo (LEMC) simulation technique and is described in Sec.(2) in detail. The method, called NP+LEMC, was applied for var- ious problems in the last couple of years. These prob- lems include particle transport through model membranes [1, 2], ion channels [3–5], and nanopores [6–8]. The ba- sic goal is to have a computationally efficient simulation technique using reduced models of steady-state systems, where particles (ions being the focus) diffuse as a result of a maintained driving force that is a concentration gra- dient and/or an electrical potential gradient (that is, an electrochemical potential gradient). We tend to regard these kind of systems as nanode- ∗Correspondence: fertig.david92@gmail.com vices that yield a macroscopic output signal (electrical current, for example) as a response to an input signal (voltage, for example). Reduced models have the advan- tage of including those degrees of freedom of the many- particle system that are essential in reproducing device behavior, namely, the relationship of the output and input signals [6]. We present two case studies here to demonstrate both the effectiveness (and also possible source of errors) of the method and the power of reduced models. In the Ryanodine Receptor (RyR) calcium release ion channel, we needed to model only those amino acids of the chan- nel protein that are inside or near the selectivity filter and hang into the ionic pathway. In the bipolar nanopore, we needed to reproduce the variation of the electric field along the pore axis properly. The most important reduc- tion in the degrees of freedom, however, is that we model water as a dielectric continuum. The effect of using im- plicit solvent instead of an explicit solvent model was dis- cussed in our previous work [6] through comparisons to molecular dynamics (MD) simulations. Our computational method is not the only one that is able to determine ionic current through biological and synthetic nanopores using reduced models. While the Brownian Dynamics (BD) simulation method [9–14] is the obvious candidate to simulate this problem, it has disadvantages from the point of view of sampling the flux of ions. The Poisson-Nernst-Planck (PNP) theory [11, 15–18] also uses the NP equation to compute flux, but uses a mean-field approach, the Poisson-Boltzmann 74 FERTIG, MÁDAI, VALISKÓ, AND BODA (PB) theory, to describe the statistical mechanical rela- tionship between the concentration and electrochemical potential profiles. To use something more state-of-the-art to provide this relationship beyond the PB level is the essence of our approach. We suggest the LEMC technique, which is a particle simulation method based on a three-dimensional model producing all the ionic correlations missing from PB. Gillespie et al. proposed a different technique, in which a density functional theory (DFT) was coupled to the NP equation [19, 20] and used to method (NP+DFT) to describe the behavior of the same ion channel that we study here, the RyR ion channel [21–23]. DFT is a con- tinuum theory, but a very sophisticated one, where ionic correlations are included through a series expansion of the free energy functional and the finite size of ions is taken into account through the fundamental measure the- ory [24]. Its accuracy was demonstrated by computing electrical double layers in extreme conditions and com- paring to Monte Carlo (MC) simulations [25, 26]. The disadvantage of this method is that it can be used effi- ciently only in one dimension. This reduction in dimen- sionality, however, is often feasible provided that we re- duce the model in an intelligent way. 2. Modeling and methodology A good model should be able to explain complex phe- nomena in a simplified way while it stays in agreement with experimental data. The nanopore systems (both natural and synthetic) studied in this work have some common features. The primitive model of electrolytes, that represents ions as charged hard spheres, is used. The interaction potential between two ions is defined by Coulomb’s law in a dielectric material: uij(r) =    ∞ if r < Ri +Rj 1 4πε0ε qiqj r if r ≥ Ri +Rj (1) where Ri and Rj are the radii, qi and qj are the charges of the different ionic species, ε0 is the permittivity of vac- uum, ε = 78.5 is the relative permittivity of the solvent, and r is the distance between two ions. Solvent (water, in this work) is not modeled explicitly; it is rather rep- resented by two response functions. Energetically, sol- vent acts as a dielectric background that screens the elec- tric field of ions by just dividing by ε in the Coulomb- potential (Eq.(1)). Dynamically, solvent molecules col- lide with ions and impede their diffusion; we take that effect into account by a diffusion coefficient function, Di(r). The ions of the electrolyte diffuse through a pore be- tween two bulk containers (see Fig.1). The pore pen- etrates a membrane that separates the two containers. The pore and the membrane are also modeled in a re- duced way; they are confined by hard walls with which the hard sphere ion cannot overlap. Other details of the Figure 1: Simulation cell used in the application of the NP+LEMC method. The model is rotationally symmetric: the three-dimensional model is obtained by rotating the figure about the z-axis. The domain confined by the blue line is the non-equilibrium transport region (the solution domain of the NP+LEMC system) that includes the pore and the two access regions. The electrochemical poten- tial in this region is not constant, but changes from one constant value to another (the values in the two bulk re- gions) in a monotonic manner to provide the driving force of ionic transport. Parameters R and H characterize sys- tem size. pores (amino acid side chains and surface charges) will be presented in Sec.(3) for the respective ion channel and nanopore systems. Ion transport is steady-state; concen- trations and electrical potentials, therefore, are kept fixed at the boundaries of the two baths on the two sides of the membrane (blue line in Fig.1). We assume that ion transport is described by drift- diffusion. Our procedure is a hybrid method that sepa- rates the configuration degrees of freedom of ions (ion positions) from the kinetic degrees of freedom (ion ve- locities). In BD simulations they are treated together and ionic current is measured by counting ions that pass through the pore. The BD method has the disadvantage of weak sampling of flux especially at low concentrations (as in our ion channel example) and in cases, when cur- rent is limited by ionic depletion zones (as in our bipolar nanopore example). Separation of the two parts of the phase space is a usual practice in equilibrium statistical mechanics, where the kinetic part is described by the ideal gas equations, while the excess quantities that produce all the peculiari- ties beyond ideal-gas behavior are from the potential en- ergy, that, in turn, depends only on the configuration coor- dinates [27–30]. Out of equilibrium, however, this sepa- ration is not so obvious. Classical mechanics is still valid, so particles’ motion can be described by Newton’s equa- tions of motion (MD simulations can be used). In the con- Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 75 tinuum solvent framework, Langevin’s equation is used in BD simulations. The meaning of various thermodynamic quantities, however, loses the solid ground that it has in equilib- rium. There is no well-established non-equilibrium sta- tistical mechanics finding its way into textbooks despite the fact that many authors made considerable advances in this field [31–36]. The most problematic quantities are those containing entropic effects, such as the chem- ical potential, µi. In the case of charged particles, we use the term electrochemical potential, but we talk about the same thing. This is a crucial quantity, because its ho- mogeneity defines thermodynamic equilibrium (let us as- sume that the other two intensive parameters, temperature and pressure, are constant). If µi is not constant, species i will diffuse until it becomes constant according to the second law of thermodynamics. The gradient of µi(r), therefore, is the driving force of particle diffusion. The empirical transport equation that is most widely used to describe transport of charged particles (electrod- iffusion) is the NP equation, ji(r) = − 1 kT Di(r)ci(r)∇µi(r), (2) where Di(r) is the diffusion coefficient profile of ionic species i, ci(r) is the concentration profile, µi(r) is the electrochemical potential profile, ji(r) is the particle flux density, k is Boltzmann’s constant and T = 298.15K is the temperature. The resulting ji(r) must satisfy the continuity equation: ∇ · ji(r) = 0. (3) The NP equation separates the configuration and kinetic parts of the phase space in a way that it provides flux as a function of three profiles that, in turn, depend only on configuration coordinates. Thus, we reduced our statisti- cal mechanical task to averaging over states in the con- figuration space just as we did in the case of equilibrium statistical mechanics. If we treat Di(r) as a parameter, the task is reduced to finding the proper relation between ci(r) and µi(r). Defining local concentration in a non- equilibrium situation is the same as in equilibrium: we compute the average number of particles in a small vol- ume element and divide it with the volume of that element (though we will also compute concentration in a different way on the basis of the Potential Distribution Theorem [37], see later). Eletrochemical potential, on the other hand, is more problematic. The trick in this case is that we assume lo- cal equilibrium (LE). If we divide the simulation cell into subvolumes Bα, we can characterize this subvol- ume with the electrochemical potential, µα i , and assume that this value is constant in Bα. We designate the µα i and cαi values to the mass centers of the Bα volume el- ements. The assumption of LE is also present in other methods (though not necessarily stated explicitly), where the ci(r) vs. µi(r) relationship is considered, such as in the PNP theory or in the NP+DFT method of Gillespie et al. [19, 20]. The main proposal of our technique was to use an MC method for this purpose. We assumed that the sub- volumes are open systems and they are in LE with a con- stant volume (V α), temperature (T ), and electrochemical potential (µα i ). Thus, we suggested to use Grand Canoni- cal Monte Carlo (GCMC) simulations, where particle in- sertion/deletions are attempted in the simulation cell. The acceptance probability depends on which subvolume we try to insert to (or where is the particle that we try to delete): min(1; pαi,χ(r)), where pαi,χ(r) = Nα i !(V α)χ (Nα i + χ)! exp ( − ∆U(r)− χµα i kT ) . (4) Here, Nα i is the number of ions of type i in subvolume Bα before insertion/deletion, ∆U(r) is the change of the system’s potential energy during particle insertion to po- sition r (or deletion from there), χ = 1 for insertion, and χ = −1 for deletion. The difference between this method and equilibrium GCMC is that µi is space-dependent and the acceptance criterion is referred to a given subvolume instead of the whole simulation cell. The effect of the surrounding of subvolume Bα, however, is taken into ac- count in the simulation through the energy change that includes all the interactions from other subvolumes, not only interactions between ions in Bα. Although it is tempting to view the subvolume Bα as a distinct thermodynamic system with its own ensemble of states and to include the effect of other subvolumes as an external constraint, this is not the case. The ensemble of states belongs to the whole system, because ion config- urations in subvolume Bα should be collected for every possible ion configurations of all the other subvolumes. Therefore, the independent variables of this ensemble are T and {V α, µα i }, where α and i run over the volume el- ements and particle species, respectively. For compari- son, the variables in global equilibrium are T , V , and µi, where V is the total volume and µi does not depend on space. We solve the NP+LEMC system in an iterative way. The electrochemical potential is adjusted until conserva- tion of mass (∇ · ji(r) = 0) is satisfied. The procedure can be summarized as µα i [n] LEMC −−−−→ cαi [n] NP −−→ jαi [n] ∇·j=0 −−−−→ µα i [n+ 1]. (5) The electrochemical potentials for the next iteration, µα i [n+ 1], are computed from the results of the previous iteration, cαi [n], on the basis of the divergence-theorem (also known as Gauss-Ostrogradsky’s theorem). The con- tinuity equation is converted to a surface integral: 0 = ∫ Bα ∇ · ji(r) dV = ∮ Sα ji(r) · n(r) da, (6) where volume Bα is bounded by surface Sα and n(r) denotes the normal vector pointing outward at position 45(1) pp. 73-84 (2017) 76 FERTIG, MÁDAI, VALISKÓ, AND BODA r of the surface. Every Sα surface is divided into Sαβ elements. Along these elements, Bα and Bβ are adjacent cells. It is assumed that the concentration, the gradient of the electrochemical potential, the flux density and the diffusion coefficient are constant on a surface Sαβ . They are denoted by hat: ĉαβi , ∇µ̂αβ i , ĵαβi , and D̂αβ i . The ĉαβi values are obtained from the values cαi and cβi via linear interpolation. The ∇µ̂αβ i values are also obtained from µα i and µβ i assuming linearity. Thus the integral in Eq.(6) for a given surface Sα is replaced by a sum over the surface elements that consti- tute Sα: 0 = ∑ β,Sαβ∈Sα ĵ αβ i · nαβaαβ , (7) where aαβ is the area of surface element Sαβ and nαβ is the outward normal vector in the center of Sαβ . The iteration procedure is described by the following steps: 1. An appropriately chosen initial set of electrochem- ical potentials is chosen (µα i [1]; in general: µα i [n], where [n] denotes the nth iteration). 2. Using these µα i [n] parameters as inputs, LEMC sim- ulations are performed. The resulting concentrations are denoted by cαi [n]. 3. The flux computed from the {µα i [n], c α i [n]} pair usu- ally does not satisfy Eq.(7). The next set of electro- chemical potential is calculated by assuming that the {µα i [n + 1], cαi [n]} pair does satisfy Eq.(7). If we write the value of ĵαβi as given by the NP equation into Eq.(7), we obtain 0 = ∑ β D̂αβ i ĉαβi [n]∇µ̂αβ i [n+ 1] · nαβ aαβ , (8) where β runs over all the surface elements Sαβ that constitute Sα. Eq.(8) is a system of linear equations; the unknown variables are denoted by µα,CAL i [n + 1], where CAL refers to the fact that these values come from calculations by solving Eq.(8). 4. To achive faster and more robust convergence in the case of large driving forces, the electrochemical po- tential used in the (n+ 1)th iteration is mixed from the values calculated in the (n+ 1)th iteration from Eq.(8) and the values mixed in the nth iteration: µα,MIX i [n+ 1] = biµ α,CAL i [n+ 1] + (1− bi)µ α,MIX i [n], (9) where bi is a mixing parameter that determines the ratio of mixing. If the parameter is close to 1, faster iteration can be achieved, however, it may result in the system fluctuating between local minima. Smaller bi values prevent this fluctuation at the price of making the convergence slower. 5. The input of the (n + 1)th LEMC simulation is µα,MIX i [n+ 1]. The system of linear equations contains the boundary conditions for the electrochemical potential and concen- tration via the boundary elements that are at the system’s boundaries (blue line in Fig.1). If the Sαβ face is on the system’s boundary, the values of ĉαβi and µ̂αβ i are those prescribed on that face. The electrochemical potentials on the left (L) and right (R) boundaries are computed from µL i = kT ln cLi + µEX,L i + qiΦ L (10) and µR i = kT ln cRi + µEX,R i + qiΦ R, (11) where µEX,L i and µEX,R i are the excess chemical poten- tials in the absence of an external field (determined by the Adaptive GCMC method of Malasics et al. [38]), while qiΦ L and qiΦ R are the interactions with the applied elec- trical potentials. Prescribing ΦL and ΦR on the system’s boundary means that we use an electrostatic Dirichlet boundary condition. Voltage is defined as U = ΦL−ΦR. Ultimately, we have boundary conditions for ion concen- trations on the two sides of the membrane and for the voltage (ground is on the right in this study). The energy change ∆U contains not only the interac- tions between particles (and interactions of particles with the pore), but also the interaction with an external elec- trical potential, Φappl(r). This applied potential is calcu- lated by solving Laplace’s equation ∇2Φappl(r) = 0 (12) for the system (inside the blue line) with the prescribed Dirichlet boundary condition on the system’s boundaries (ΦL and ΦR on the left and right blue line, respectively). In the portion of the blue line inside the membrane, we interpolate between ΦL and ΦR. We solve this equation with the Induced Charge Computation method [39, 40]. In the present geometry, the elementary cells are ∆z× ∆r rectangles in the (z, r) plane. In three-dimensional space, these correspond to concentric rings with the z- axis in their centers. The concentration in an elementary cell Bα can be computed from dividing the average number of ions in the cell with the volume of the cell. This route is disad- vantageous if ion concentrations are small. An alternative method is based on the Potential Distribution Theorem [37] that practically corresponds to the Widom particle insertion method [41, 42]. The total excess chemical po- tential (containing also the interaction with the applied field) is µα i − kT ln cαi and can be computed as exp [− (µα i /kT − ln cαi )] = 〈exp [−∆Uα i /kT ]〉 , (13) where ∆Uα i is the energy change associated with the in- sertion of a test particle of species i into a randomly cho- sen position in the elementary cell α. This is the same en- ergy computed in the particle insertion step of the LEMC technique, therefore, it does not require additional com- putational cost. The concentration can be calculated from Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 77 Eq.(13) because µα i is known (that is the input of the sim- ulation) and the ensemble average on the right hand side can be obtained from the LEMC insertions. This route provides a more accurate value for cαi , because this sam- pling is not discrete as just counting particles, but rather continuous, because we always gain information from the simulation at every ion insertion via the energy ∆Uα i . The route of counting particles, however, might be bet- ter, when concentrations are very large in the volume el- ement. In the case of long enough simulations, the two methods give identical results (within a statistical error). Because of this statistical error, the iteration does not converge to an exact value of the ionic flux density, but it fluctuates around a limiting value. The final solu- tion is obtained from a running average over iterations. Longer LEMC simulations and more iterations result in a more reliable outcome. The net flux of the diffusing ions through the cross section of the pore is calculated via av- eraging: Ji(z) = 2π ∫ R(z) 0 rji(z, r) · nz dr (14) where R(z) denotes the radius of the pore at coordinate z (|z| < Hmemb/2, where Hmemb is the thickness of the membrane) and nz is the unit vector along the z-axis. The electrical current of an ion is then Ii = qiJi, (15) and the total current is I = ∑ i qiJi. (16) 3. Results and discussion 3.1 Calcium release channel There are two important classes of calcium channels that have been our focus. The L-type calcium channel is found in excitable cell membranes [43] such as those of nerve cells and muscle cells (both cardiac and skeletal). These are strongly Ca2+ selective channels, meaning that the presence of only a µM Ca2+ in the bulk decreases the Na+ current to half of its value compared to its value in the absence of Ca2+. This is called micromolar Ca2+ affinity (or Ca2+ block) because a very small amount of Ca2+ is enough to block the channel [44]. The price of high Ca2+ selectivity is low Ca2+ current through this channel. The Ca2+ ions flow into the muscle cell when the L-type calcium channels open in response to an electrical signal (action potential) in a small quan- tity that is not sufficient to initiate muscle contraction. The solution of Nature for this problem is an amplifica- tion mechanism. A large amount of Ca2+ ions are stored in organelles, the sarcoplasmic reticulum (SR), that are situated inside the muscle cells. There are calcium release channels (also known as RyR calcium channels) embedded in the mem- brane of the SR that are opened by the Ca2+ ions pro- vided by the L-type channel (in cardiac muscle) or via a direct link between the two types of calcium channels (in skeletal muscle). The RyR channels provide the large Ca2+ flux that is necessary for muscle contraction by binding to tropomyosins on the actin filaments and al- lowing myosin to climb up the filament. These channels are wider (also, less Ca2+ selective) than the L-type channels. A large amount of experimen- tal data is available for this channel [45–47]. On the basis of this database, Gillespie et al. developed a model for the RyR channel [21,22]. Although the model was a one- dimensional reduction, studied with a one-dimensional NP+DFT, they were able to reproduce hundreds of ex- perimental current-voltage curves. Later, when the NP+LEMC method became avail- able, we developed [4] a three-dimensional version of the model of Gillespie et al. The model is similar to that used by Boda et al. [3, 4, 40, 48–54] for the L-type cal- cium channel inspired by the “space-charge competition mechanism” proposed by Nonner et al. [55]. From the point of view of ion selectivity and permeation, the im- portant part of the channel is the selectivity filter, the bot- tleneck of the pore. This region is lined by four P-loops that are parts of the four membrane-spanning subunits. These loops contribute four glutamic acids (in the L-type calcium channel) to the filter-lining region (the EEEE lo- cus). The COO− groups of the carboxyl side chains are thought to reach into the pore lumen and interact with the passing ions [44]. In the case of the RyR channel, they are four aspartates (D4899, see Fig.2). Gillespie et al. [21–23, 46, 47] identified additional E and D amino acids in the two vestibules on the cytosolic and luminal sides (Fig.2). It has been shown [3, 4, 40, 48–55] that the strong Ca2+ selectivity can be reproduced by a reduced model where eight half-charged oxygen atoms (O1/2−) of the four COO− groups force a strict competition between Ca2+ and Na+ ions for space in the crowded filter. The exact positions of the O1/2− ions are not known and even irrelevant from the point of view of the selectivity of the model. What matters is that they attract the cations into the filter with their charge, and, at the same time, they ex- ert a hard sphere exclusion and make it difficult for the cations to find space in the filter. Therefore, it proved to be sufficient to model the structural groups of the filter as mobile O1/2− ions that are confined to the filter re- gion. The profiles shown in Fig.2 are results of the sim- ulations and show how the eight O1/2− ions distribute in a given segment of the channel. There are 32 such O1/2− ions, plus point charges placed on a ring on the luminal side. The geometrical features (filter radius is 5.5 Å, filter length is 15 Å, and radii of left and right vestibules are 22 and 9 Å, respectively) are also designed on the basis of Gillespie’s model that, in turn, are based on fitting to experimental data. This model was designed and used unchanged in ev- 45(1) pp. 73-84 (2017) 78 FERTIG, MÁDAI, VALISKÓ, AND BODA Figure 2: The three-dimensional model [4] of the RyR calcium channel based on the model of Gillespie et al. [21,22]. Each curve shows the distribution of eight O1/2− ions that are confined to the given region but otherwise free to move there. The charges of the E4902 amino acids are represented as point charges (−1/2e) at fixed positions on a ring. The thickness of the membrane is Hmemb = 46 Å. The dielectric constant is ε = 78.5. ery calculation. What remains to be specified are the dif- fusion coefficient profiles of the various ionic species. These profiles depend only on z, Di(z); we assume no variation in the radial dimension. The values outside the channel and in the selectivity filter, DB i and DF i , respec- tively, are constant. In the vestibules at the two entrances of the channel, the diffusion constant profiles are interpo- lated between the DB i and DF i values in a way that their value changes in proportion with the cross section of the channel. The bulk vales, DB i , are taken from experiments; they are 1.334 ·10−9, 7.92 ·10−10, and 2.032 ·10−9 m2s−1 for Na+, Ca2+, and Cl−, respectively. The values inside the cylindrical selectivity filter are fitted to two experimental data points: 250 mM symmetric NaCl at 100 mV for Na+, while added 10 mM luminal CaCl2 at -100 mV for Ca2+. The resulting values (1.27·10−10 and 1.27·10−11 for Na+ and Ca2+, respectively) are fixed and never changed. The Cl− value is irrelevant, because Cl− ions do not carry significant current. These values have been obtained for a moderate system size (H = 54 Å, R = 48 Å). Currents, and, therefore, diffusion coefficients fitted to currents can depend on system size (see later). Here, we show results for a NaCl-CaCl2 mixture, where there is 250 mM Na+ on both sides of the mem- brane, while there is 4 µM Ca2+ on the cytosolic (left) and 50 mM Ca2+ on the luminal (right) side. The volt- age is changed from -150 mV to 150 mV with the ground on the right. The current vs. voltage curves are shown Figure 3: Currents vs. voltage curves as obtained from ex- periment [21], the model of Gillespie et al. [21, 22] ob- tained from DFT coupled to the NP equation, and from the NP+LEMC method. In the case of the NP+LEMC method, the currents carried by Na+ and Ca2+ are also shown. in Fig.3. Total currents are shown in black as obtained from experiments (× symbols), Gillespie’s NP+DFT cal- culations (dashed line), and our NP+LEMC calculations (solid line with filled triangles). Agreement with experiment is very good using both models. The slope of the I − U curve is larger for posi- tive voltages than for negative voltages. More details and understanding can be gained from the current curves for the separate ions, Na+ and Ca2+ (blue and red curves, respectively). At positive voltages, the Ca2+ current is practically zero, so the total current is carried by Na+ ions. At negative voltages, on the other hand, both Na+ and Ca2+ ions contribute to the current. The explanation of this behavior can be drawn from concentration and electrochemical potential profiles shown in Fig.4. Let us discuss the Na+-profiles first (blue curves). The driving force for the Na+ ions is the voltage because Na+ concentration is the same on the two sides. The difference between Na+ currents at -100 and 100 mV voltages, therefore, is the result of the different competi- tion between Na+ and Ca2+ ions at the two voltages. To understand the difference in this competition, we need to understand the behavior of Ca2+ ions. At 100 mV, the driving force for Ca2+ ions is small because the concentration difference balances the electri- cal potential difference (see the Ca2+ electrochemical po- tential profile in the right panel of Fig.4B). This results in a small Ca2+ current. Ca2+ concentration is small in the left bulk, therefore, a Ca2+ depletion zone is formed in the left vestibule and in the selectivity filter of the channel (beware the logarithmic scale of the concentration axis). That depletion zone results in a decreased concentration of Ca2+ ions compared to Na+ ions inside the channel. In the case of -100 mV, on the other hand, the Ca2+ depletion zone in the left vestibule is absent because Ca2+ ions arrive from the right and “fill up” the left Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 79 (A) (B) Figure 4: (A) Concentration and (B) electrochemical po- tential profiles of Na+ and Ca2+ for two opposing volt- ages: -100 mV (left panels) and 100 mV (right pan- els). The striped parts indicate the cytosolic and luminal vestibules, while the gray area is the selectivity filter. vestibule. There is a more balanced competition between Ca2+ and Na+ ions in the channel and Ca2+ ions can use the free-energetic advantage that they have over Na+ [22, 53]. This means that Ca+ ions are favored by the crowded selectivity filter, because they provide twice the charge (compared to Na+ ions) to balance the charge of O1/2− ions while occupying about the same space (their diameters are similar: 1.98 and 1.9 Å for Ca2+ and Na+, respectively). The finite size of the ions plays a crucial role in the selectivity mechanism. This kind of selectivity could not be produced with the PNP theory. Because Ca2+ concentration is large in the left vestibule of the channel, it drops quickly to the 4 µM value at the left boundary of the solution domain. This in- troduces a severe system-size dependence into the calcu- lations in this case that is analyzed in Fig.5. Small Ca2+ concentration on the left hand side corresponds to a large Debye-length in the double layer at left hand side of the membrane (in the left access region). That double layer should fit into the simulation cell (as it does in the case of a larger cell, see black curves in Fig.5). If the cell is too small (see red curves in Fig.5), there is not enough space for the Ca2+ concentration to reach the 4 µM limiting value at the left boundary of the cell. The electrochemical potential cannot reach its limiting value either (see the bottom panel). The NP+LEMC calcula- tion, however, provides a solution in this case too, be- cause the layer near the left outer boundary of the cell takes care of the missing access region in an averaged manner. The concentration has a small value in the layer, Figure 5: Concentration (top panel) and electrochemical potential (bottom panel) profiles of Ca2+ for -100 mV. The results of simulations for two system sizes are shown: H = 54 Å (red) and H = 180 Å (black). while the electrochemical potential profile drops abruptly (large driving force). An appropriate value of the ci∇µi product, therefore, is provided by the self-consistent so- lution so that the continuity equation is satisfied. This solution, however, is approximate. A large por- tion of the access region with considerable resistance is taken into account in an averaged way by the layer near the system edge. This introduces an error into the value of the Ca2+ current that is indicated in Fig.5 for both cases. The good behavior of the current-voltage curves compared to experiments and DFT calculations is due to the fit of the Ca2+ diffusion coefficient, DF Ca2+ , inside the filter. That value was fitted for negative voltage at the given system size balancing the system-size error. This result points out the importance of system size in the case of small ionic concentrations, but it also shows the role of the diffusion coefficient profile in NP+LEMC calculations. In confined geometries, the Di(r) profile, although it has a strong relation to the mobility of ions, is primarily an adjustable parameter that is fitted to experi- ments (as in the present case) or to MD results [6]. The value of DF i takes into account interactions that are ab- sent in the reduced model or accounts for resistances of regions that are absent in the model. In the case of the RyR channel, for example, the dif- fusion coefficient profile includes effects of the parts of the ion channel not included in the model: the real RyR channel is much larger than the 46 Å portion modeled here. That region also tunes the total current, but it is not selective. The selectivity filter and its close neigh- 45(1) pp. 73-84 (2017) 80 FERTIG, MÁDAI, VALISKÓ, AND BODA Figure 6: Bipolar nanopore geometries: Cyl: cylindrical; SC: single conical; DC: double conical. The pore is pos- itively charged on the left hand side (−30 < z < 0 Å), while negatively charged on the right hand side (0 < z < 30 Å) keeping the surface charge fixed (σ = ±0.1 C/m2). Minimal and maximal pore radii are 10 and 20 Å, respec- tively. borhood modeled here determines ion selectivity and is able to reproduce complex behavior such as anomalous mole fraction experiments discussed previously by Gille- spie [21–23] and Boda [4]. 3.2 Rectifying bipolar nanopores Ion channels are natural nanopores with stable well- defined structures that are very narrow at their selectivity filters (often below 1 nm in radius) to make them suit- ably selective. The disadvantages, however, are consider- able. The structures are often unknown. They are difficult to handle experimentally. Their manipulation is cumber- some with point mutations. Synthetic nanopores, therefore, quickly gained atten- tion due to the fact that they have special properties com- pared to those of micropores. These special properties arise because the screening length of the electrolyte is comparable to the radius of the nanopore. This fact gives nanopores properties that resemble those of ion channels. One advantage of synthetic nanopores is that are rela- tively easy to fabricate [56–63]. They are either solid state nanopores using ion-beam or electron-beam sculpting in, for example, silicon compound membranes, or they are track etched into polymer membranes. Two basic proper- ties of such nanopores are their geometry (shape) and the pattern of surface charge on the pore wall. Figure 7: Current-voltage relations for the three nanopore geometries (Cyl, SC, and DC from bottom to top). Cur- rents carried by Na+ (solid blue), Cl− (dashed red), and their sum (dot-dashed black) are shown as a function of voltage. The insets magnify the results for negative volt- ages. In our previous works, we studied the bipolar nanopore, where the surface charge is positive on the left hand side of the pore, while it is negative on the right hand side [6, 7]. These nanopores rectify ionic current, mean- ing that they let a much larger amount of ions through at a given sign of voltage than at the opposite sign. In those papers, we used a cylindrical geometry for the nanopore and focused on the effect of the charge pattern, pore ra- dius, and concentration. Here, we discuss the effect of pore geometry on the rectification properties of a bipolar nanopore. We per- formed NP+LEMC calculations for three different ge- ometries, the cylindrical (Cyl), single conical (SC), and double conical (DC) shown in Fig.6. Simulations for dif- ferent voltages have been performed from -150 mV to 150 mV. The concentration of the electrolyte was 0.1 M on both sides. The electrolyte is a 1:1 system (let us call it NaCl), but the diameters of the cations and anions are the same in order to avoid effects from ion size asymmetry. For the same reason, the same diffusion coefficients were used for the two ions. This makes it possible to focus on the balance of charge and geometrical asymmetries. If the nanopore’s shape is symmetric, for example, the cations and anions carry the same amount of current (Cyl and DC). Figure 7 shows the current-voltage relations. Rectifi- cation is observed in all three geometries: current is much larger at positive than at negative voltages. The rectifi- Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 81 Figure 8: Rectification (defined as the absolute value of the current ratio in the ON and OFF states) as a function of the absolute value of the voltage for the three nanopore geometries. cation is defined as the ratio of currents at positive (ON state) and negative voltages (OFF state). The results are shown in Fig.8. The basic explanation of rectification can be depicted from the concentration profiles (Fig.9). There are deple- tion zones for cations in the positively charged half region (left), while there are depletion zones for anions in the negatively charged half region (right). In the OFF state (negative voltage) these depletion zones are deeper than in the ON state (positive voltage). More detailed expla- nation have been given in our previous works [6,7]. Here we focus our discussion to the effect of pore shape. Total currents are larger in the SC and DC geometries due to their wider entrances. Interestingly, total currents are the same in the SC and DC geometries despite the quite different pore shapes. Whether this is a coincidence or it has a deeper explanation requires further investiga- tion. What is different in the SC and DC geometries is the partitioning of the total current between Na+ and Cl−. While Na+ and Cl− currents are the same in the DC geometry due to the symmetric shape, they are differ- ent in the SC geometry. The current carried by Na+ ions is smaller than the current carried by Cl− ions (middle panel of Fig.7). This is reflected in concentration profiles (see middle panel of Fig.9). The explanation is that the tip of the SC nanopore (left entrance) is positively charged so the depletion zones of Na+ ions there is deeper than the depletion zone of Cl− ions on the other side where the pore is wider. Deeper depletion zones result in smaller currents. Therefore, Na+ current is smaller in the whole voltage range, but more so in the OFF state at negative voltages. Rectification is largest in the Cyl geometry because the depletion zones are the deepest in that geometry for both ions. The deepest points of the depletion zones are formed around |z| = 10 Å. There are different effects that form the concentration profiles inside the pore. In the left region, for example, (1) the positive surface charge on the pore wall repels Na+ ions, (2) the negative surface charge Figure 9: Concentration profiles of Na+ (blue) and Cl− (red) in the ON (solid lines) and OFF (dashed lines) states for the three nanopore geometries (Cyl, SC, and DC from bottom to top). on the right hand side attracts them, (3) the 0.1 M bulk region acts as a source for the ions, (4) applied potential in the OFF state drives Na+ ion to the right, where they have a peak, and (5) Cl− ions (that have a peak on the left) attract them. The balance of all these effects forms the deep depletion zone of Na+ at z ≈ −10 Å in the OFF state. Because the Cyl geometry has the smallest radius at the |z| ≈ 10 Å positions, this geometry provides the deepest depletion zone as a result of the dominant effect from the above list, the effect of repelling surface charge. 4. Summary We presented the NP+LEMC method that is a hybrid technique, harvesting the advantageous properties of both the NP transport equation (fast calculation of flux) and LEMC particle simulation method (correct calculation of ionic correlations). We applied the method to compute ionic currents through reduced models of an ionic chan- nel and a bipolar nanopore. Reduced models have the ad- vantage of fast calculation (LEMC would not be feasible for an explicit-water model) and intellectual focus. With these models, we can concentrate on those properties of the system that are important to reproduce its behavior as a device. For the RyR calcium channel, for example, the re- duced model is able to reproduce complex selectivity behavior in agreement with experiments by modeling only the “important” amino acids in a reduced way 45(1) pp. 73-84 (2017) 82 FERTIG, MÁDAI, VALISKÓ, AND BODA [21,22]. For the bipolar nanopore, the reduced model us- ing implicit water is able to reproduce MD results for an explicit-water model [6]. With further methodological and model development, we intend to simulate nanode- vices as close to their real size as possible. Acknowledgement The financial support of the National Research, Develop- ment and Innovation Office (NKFIH K124353) is grate- fully acknowledged. The present article was published in the frame of the Project No. GINOP-2.3.2-15-2016- 00053 and supported by the UNKP-17-4 (to M.V.) and UNKP-17-1 (to E.M.), the New National Excellence Pro- gram of the Ministry of Human Capacities. We thank Dirk Gillespie for helpful discussions. REFERENCES [1] Boda, D., Gillespie, D.: Steady state electrodif- fusion from the Nernst-Planck equation coupled to Local Equilibrium Monte Carlo simulations, J. Chem. Theor. Comput., 2012 8(3), 824–829, DOI: 10.1021/ct2007988 [2] Ható, Z., Boda, D., Kristóf, T.: Simulation of steady-state diffusion: Driving force ensured by Dual Control Volumes or Local Equilibrium Monte Carlo, J. Chem. Phys., 2012 137(5), 054109, DOI: 10.1063/1.4739255 [3] Boda, D., Kovács, R., Gillespie, D., Kristóf, T.: Se- lective transport through a model calcium channel studied by Local Equilibrium Monte Carlo simula- tions coupled to the Nernst-Planck equation, J. Mol. Liq., 2014 189, 100, DOI: 10.1016/j.molliq.2013.03.015 [4] Boda, D.: Monte Carlo Simulation of Electrolyte Solutions in Biology: In and Out of Equilib- rium, Annual Reports in Computational Chemistry, 2014, vol. 10, chap. 5 (Elsevier), 127–163, DOI: 10.1016/b978-0-444-63378-1.00005-7 [5] Ható, Z., Boda, D., Gillepie, D., Vrabec, J., Rutkai, G., Kristóf, T.: Simulation study of a rectifying bipolar ion channel: detailed model versus reduced model, Cond. Matt. Phys., 2016 19(1), 13802, DOI: 10.5488/cmp.19.13802 [6] Ható, Z., Valiskó, M., Kristóf, T., Gillespie, D., Boda, D.: Multiscale modeling of a rectifying bipolar nanopore: explicit-water versus implicit- water simulations, Phys. Chem. Chem. Phys., 2017 17(27), 17816–17826, DOI: 10.1039/C7CP01819C [7] Matejczyk, B., Valiskó, M., Wolfram, M.T., Pietschmann, J.F., Boda, D.: Multiscale modeling of a rectifying bipolar nanopore: Comparing Poisson- Nernst-Planck to Monte Carlo, J. Chem. Phys., 2017 146(12), 124125, DOI: 10.1063/1.4978942 [8] Mádai, E., Valiskó, M., Dallos, A., Boda, D.: Sim- ulation of a model nanopore sensor: Ion competi- tion underlies device behavior, J. Chem. Phys., 2017 147(24), 244702, DOI: 10.1063/1.5007654 [9] van Gunsteren, W., Berendsen, H.: Algorithms for Brownian dynamics, Mol. Phys., 1982 45(3), 637– 647, DOI: 10.1080/00268978200100491 [10] Turq, P., Lantelme, F., Friedman, H.L.: Brown- ian dynamics: Its application to ionic solutions, J. Chem. Phys., 1977 66(7), 3039–3044, DOI: 10.1063/1.434317 [11] Corry, B., Kuyucak, S., Chung, S.H.: Tests of con- tinuum theories as models of ion channels. II. Poisson-Nernst-Planck theory versus Brownian dy- namics, Biophys. J., 2000 78(5), 2364–2381, DOI: 10.1016/s0006-3495(00)76781-6 [12] Chung, S.H., Allen, T.W., Hoyles, M., Kuyucak, S.: Permeation Of Ions Across The Potassium Chan- nel: Brownian Dynamics Studies, Biophys. J., 1999 77(5), 2517–2533, DOI: 10.1016/s0006-3495/99/77087-6 [13] Chung, S.H., Hoyles, M., Allen, T., Kuyucak, S.: Study Of Ionic Currents Across A Model Mem- brane Channel Using Brownian Dynamics, Bio- phys. J., 1998 75(2), 793–809, DOI: 10.1016/s0006- 3495(98)77569-1 [14] Berti, C., Furini, S., Gillespie, D., Boda, D., Eisen- berg, R.S., Sangiorgi, E., Fiegna, C.: A 3-D Brow- nian Dynamics simulator for the study of ion permeation through membrane pores, J. Chem. Theor. Comput., 2014 10(8), 2911–2926, DOI: 10.1021/ct4011008 [15] Zheng, Q., Chen, D., Wei, G.W.: Second-order Poisson-Nernst-Planck solver for ion transport, J. Comp. Phys., 2011 230(13), 5239–5262, DOI: 10.1016/j.jcp.2011.03.020 [16] Schuss, Z., Nadler, B., Eisenberg, R.S.: Deriva- tion of Poisson and Nernst-Planck equations in a bath and channel from a molecular model, Phys. Rev. E, 2001 6403(3), 036116, DOI: 10.1103/phys- reve.64.036116 [17] Nonner, W., Eisenberg, B.: Ion permeation and glu- tamate residues linked by Poisson-Nernst-Planck theory in L-type calcium channels, Biophys. J., 1998 75(3), 1287–1305, DOI: 10.1016/s0006- 3495(98)74048-2 [18] Pietschmann, J.F., Wolfram, M.T., Burger, M., Trautmann, C., Nguyen, G., Pevarnik, M., Bayer, V., Siwy, Z.: Rectification properties of conically shaped nanopores: consequences of miniaturiza- tion, Phys. Chem. Chem. Phys., 2013 15(3), 16917– 16926, DOI: 10.1039/C3CP53105H [19] Gillespie, D., Nonner, W., Eisenberg, R.S.: Cou- pling Poisson-Nernst-Planck and density functional theory to calculate ion flux, J. Phys.: Cond. Matt., 2002 14(46), 12129–12145, DOI: 10.1088/0953- 8984/14/46/317 [20] Gillespie, D., Nonner, W., Eisenberg, R.S.: Density functional theory of charged, hard-sphere fluids, Phys. Rev. E, 2003 68(3), 031503, DOI: 10.1103/phys- reve.68.031503 [21] Gillespie, D., Xu, L., Wang, Y., Meissner, G.: (De)constructing the ryanodine receptor: Modeling Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 83 ion permeation and selectivity of the calcium re- lease channel, J. Phys. Chem. B, 2005 109(32), 15598–15610, DOI: 10.1021/jp052471j [22] Gillespie, D.: Energetics of Divalent Selectivity in a Calcium Channel: The Ryanodine Receptor Case Study, Biophys. J., 2008 94(4), 1169–1184, DOI: 10.1529/biophysj.107.116798 [23] Gillespie, D., Giri, J., Fill, M.: Reinterpreting the anomalous mole fraction effect: The Ryanodine re- ceptor case study, Biophys. J., 2009 97(8), 2212– 2221, DOI: 10.1016/j.bpj.2009.08.009 [24] Rosenfeld, Y.: Free-energy model for the inho- mogeneous hard-sphere fluid mixture and density- functional theory of freezing, Phys. Rev. Lett., 1989 63(9), 980–983, DOI: 10.1103/physrevlett.63.980 [25] Gillespie, D., Valiskó, M., Boda, D.: Density func- tional theory of the electrical double layer: the RFD functional, J. Phys.-Cond. Matt., 2005 17(42), 6609–6626, DOI: 10.1088/0953-8984/17/42/002 [26] Valiskó, M., Gillespie, D., Boda, D.: Selective ad- sorption of ions with different diameter and valence at highly-charged interfaces, J. Phys. Chem. C, 2007 111(43), 15575–15585, DOI: 10.1021/jp073703c [27] McQuarrie, D.A.: Statistical Mechanics (Univer- sity Science Books, Sausalito), 2000, ISBN: 9781891389153 [28] Allen, M.P., Tildesley, D.J.: Computer Simulation of Liquids, Oxford Science Publ (Clarendon Press, New York), 1989, ISBN: 9780198556459 [29] Frenkel, D., Smit, B.: Understanding molecular simulations (Academic Press, San Diego), 1996, ISBN: 9780122673702 [30] Hansen, J.P., McDonald, I.R.: Theory of Simple Liquids (Elsevier, Academic Press), 2006 ISBN: 9780080455075 [31] Prigogine, I.: Introduction to Thermodynamics of Irreversible Processes (John Wiley and Sons, New York), 1968, ISBN: 9780470699287 [32] Kreuzer, H.J.: Nonequilibrium Thermodynamics and its Statistical Foundations, Monographs on the physics and chemistry of materials (Clarendon Press, Oxford), 1983, ISBN: 9780198513759 [33] de Groot, S.R., Mazur, P.: Non-Equilibrium Ther- modynamics, Dover Books on Physics (Dover, New York), 2013, ISBN: 9780486153506 [34] Tschoegl, N.W.: Fundamentals of Equilibrium and Steady-State Thermodynamics (Elsevier, Amster- dam), 2000, ISBN: 9780080532110 [35] Attard, P.: Thermodynamics and Statistical Mechanics, Equilibrium by Entropy Maximi- sation (Academic Press, London), 2002, ISBN: 9780120663217 [36] Evans, D.J., Morriss, G.: Statistical Mechanics of Nonequilibrium Liquids (Cambridge University Press, New York), 2008, ISBN: 9781139471930 [37] Beck, T.L., Paulaitis, M.E., Pratt, L.R.: The Poten- tial Distribution Theorem and Models of Molecu- lar Solutions (Cambridge University Press, Cam- bridge), 2006, ISBN: 9781139457637 [38] Malasics, A., Boda, D.: An efficient iterative grand canonical Monte Carlo algorithm to determine in- dividual ionic chemical potentials in electrolytes, J. Chem. Phys., 2010 132(24), 244103, DOI: 10.1063/1.3443558 [39] Boda, D., Gillespie, D., Nonner, W., Henderson, D., Eisenberg, B.: Computing induced charges in inhomogeneous dielectric media: Application in a Monte Carlo simulation of complex ionic systems, Phys. Rev. E, 2004 69(4), 046702, DOI: 10.1103/phys- reve.69.046702 [40] Boda, D., Valiskó, M., Eisenberg, B., Nonner, W., Henderson, D., Gillespie, D.: The effect of protein dielectric coefficient on the ionic selectivity of a cal- cium channel, J. Chem. Phys., 2006 125(3), 034901, DOI: 10.1063/1.2212423 [41] Widom, B.: Some Topics in the Theory of Flu- ids, J. Chem. Phys., 1963 39(11), 2808–2812, DOI: 10.1063/1.1734110 [42] Widom, B.: Structure of interfaces from uniformity of the chemical potential, J. Stat. Phys., 1978 19(6), 563–574, DOI: 10.1007/bf01011768 [43] Hille, B.: Ion channels of excitable membranes (Sinauer Associates, Sunderland), 2001, ISBN: 9780878933211 [44] Sather, W.A., McCleskey, E.W.: Permeation and selectivity in calcium channels, Ann. Rev. Phys- iology, 2003 65(1), 133–159, DOI: 10.1146/an- nurev.physiol.65.092101.142345 [45] Gao, L., Balshaw, D., Xu, L., Tripathy, A., Xin, C., Meissner, G.: Evidence for a Role of the Lume- nal M3-M4 Loop in Skeletal Muscle Ca2+ Release Channel (Ryanodine Receptor) Activity and Con- ductance, 2000 79(2), 828–840, DOI: 10.1016/s0006- 3495(00)76339-9 [46] Wang, Y., Xu, L., Pasek, D.A., Gillespie, D., Meiss- ner, G.: Probing the Role of Negatively Charged Amino Acid Residues in Ion Permeation of Skeletal Muscle Ryanodine Receptor, 2005 89(1), 256–265, DOI: 10.1529/biophysj.104.056002 [47] Xu, L., Wang, Y., Gillespie, D., Meissner, G.: Two Rings of Negative Charges in the Cytosolic Vestibule of Type-1 Ryanodine Receptor Modulate Ion Fluxes, 2006 90(2), 443–453, DOI: 10.1529/bio- physj.105.072538 [48] Boda, D., Valiskó, M., Eisenberg, B., Nonner, W., Henderson, D., Gillespie, D.: Combined ef- fect of pore radius and protein dielectric coeffi- cient on the selectivity of a calcium channel, Phys. Rev. Lett., 2007 98(16), 168102, DOI: 10.1103/phys- revlett.98.168102 [49] Gillespie, D., Boda, D.: The anomalous mole frac- tion effect in calcium channels: A measure of pref- erential selectivity, Biophys. J., 2008 95(6), 2658– 2672, DOI: 10.1529/biophysj.107.127977 45(1) pp. 73-84 (2017) 84 FERTIG, MÁDAI, VALISKÓ, AND BODA [50] Boda, D., Valiskó, M., Henderson, D., Eisenberg, B., Gillespie, D., Nonner, W.: Ion selectivity in L- type calcium channels by electrostatics and hard- core repulsion, J. Gen. Physiol., 2009 133(5), 497– 509, DOI: 10.1085/jgp.200910211 [51] Malasics, M., Boda, D., Valiskó, M., Henderson, D., Gillespie, D.: Simulations of calcium channel block by trivalent ions: Gd3+ competes with per- meant ions for the selectivity filter, Biochim. et Bio- phys. Acta - Biomembranes, 2010 1798(11), 2013– 2021, DOI: 10.1016/j.bbamem.2010.08.001 [52] Rutkai, G., Boda, D., Kristóf, T.: Relating bind- ing affinity to dynamical selectivity from dynamic Monte Carlo simulations of a model calcium chan- nel, J. Phys. Chem. Lett., 2010 1(14), 2179–2184, DOI: 10.1021/jz100718n [53] Boda, D., Giri, J., Henderson, D., Eisenberg, B., Gillespie, D.: Analyzing the components of the free energy landscape in a calcium selective ion chan- nel by Widom’s particle insertion method, J. Chem. Phys., 2011 134(5), 055102, DOI: 10.1063/1.3532937 [54] Boda, D., Henderson, D., Gillespie, D.: The role of solvation in the binding selectivity of the L-type cal- cium channel, J. Chem. Phys., 2013 139(5), 055103, DOI: 10.1063/1.4817205 [55] Nonner, W., Catacuzzeno, L., Eisenberg, B.: Bind- ing and selectivity in L-type calcium channels: A mean spherical approximation, Biophys. J., 2000 79(4), 1976–1992, DOI: 10.1016/s0006-3495(00)76446-0 [56] Siwy, Z., Fulinski, A.: Fabrication of a synthetic nanopore ion pump, Phys. Rev. Lett., 2002 89(19), 198103, DOI: 10.1103/physrevlett.89.198103 [57] Siwy, Z., Apel, P., Baur, D., Dobrev, D.D., Korchev, Y.E., Neumann, R., Spohr, R., Trautmann, C., Voss, K.O.: Preparation of synthetic nanopores with trans- port properties analogous to biological channels, Surf. Sci., 2003 532, 1061–1066, DOI: 10.1016/s0039- 6028(03)00448-5 [58] Siwy, Z., Apel, P., Dobrev, D., Neumann, R., Spohr, R., Trautmann, C., Voss, K.: Ion transport through asymmetric nanopores prepared by ion track etch- ing, Nuclear Instruments & Methods In Phys. Re- search Section B-beam Interactions With Materi- als Atoms, 2003 208, 143–148, DOI: 10.1016/s0168- 583x(03)00884-x [59] Howorka, S., Siwy, Z.: Nanopores: Generation, En- gineering, and Singly-Molecule Applications, in Handbook of Single-Molecule Biophysics (P. Hin- terdorfer, A. van Oijen, eds.), Advances in Chemi- cal Physics, chap. Chapter 11 (Springer), 2009 293– 339, DOI: 10.1007/978-0-387-76497-9_11 [60] Guo, W., Tian, Y., Jiang, L.: Asymmetric Ion Transport through Ion-Channel-Mimetic Solid- State Nanopores, Acc. Chem. Res., 2013 46(12), 2834–2846, DOI: 10.1021/ar400024p [61] Gibb, T., Ayub, M.: Solid-State Nanopore Fabri- cation, in Engineered Nanopores for Bioanalytical Applications (Elsevier BV), 2013 121–140, DOI: 10.1016/b978-1-4377-3473-7.00005-4 [62] Guan, W., Li, S.X., Reed, M.A.: Voltage gated ion and molecule transport in engineered nanochan- nels: theory, fabrication and applications, Nanotech- nology, 2014 25(12), 122001, DOI: 10.1088/0957- 4484/25/12/122001 [63] Zhang, H., Tian, Y., Jiang, L.: Fundamen- tal studies and practical applications of bio- inspired smart solid-state nanopores and nanochan- nels, Nano Today, 2016 8(3), 1470–1478, DOI: 10.1016/j.nantod.2015.11.001 Hungarian Journal of Industry and Chemistry