The effect of chelating agent on hydrodesulfurization reaction of ordered mesoporous alumina supported NiMo catalysts

Di Hu ,Hui-Ping Li ,Jin-Lin Mei,Cheng-Kun Xiao,En-Hua Wang,Xi-Yue Chen,Wen-Xin Zhang,Ai-Jun Duan

State Key Laboratory of Heavy Oil Processing,China University of Petroleum-Beijing,Beijing,102249,China

Keywords:Chelating agent Hydrodesulfurization Ordered mesoporous alumina

ABSTRACT In this paper,ordered mesoporous alumina (OMA) support with the high surface area (328 m2 g-1) and the large pore volume 0.74 (cm3·g-1) was synthesized by homogeneous precipitation method.And the influence of EDTA on the physical and chemical properties of the modified catalysts was also studied.The characteristic results showed that the addition of EDTA could adjust the metal-support interaction and improved the acidity of the corresponding catalyst.Combined with the catalytic performance results,the EDTA-modified NiMoE(1.0)/OMA catalyst displays the highest DBT hydrodesulfurization conversion(97.7%).

1.Introduction

The environmental pollution brought by automobiles emission is sharply increasing.The exhaust gas generated by the combustion of sulfur compounds in gasoline and diesel is a major source of pollution to acid rain and haze weather.To produce ultra-low sulfur transportation oil,one of the key measures is to develop efficient hydrodesulfurization catalyst for sulfur removal from the fuels(Sun and Prins,2010;Escobar et al.,2018;Vatutina et al.,2016).

Traditional crystalline alumina has excellent mechanical properties,suitable acid properties,and low price,especially the most widely used γ-Al2O3(Moser et al.,2010;Santolalla-Vargas et al.,2015).But the specific surface areas of these alumina are usually less than 300 m2/g,and the pore size distributions are also relatively wide range(generally 3-15 nm).In recent years,researchers have focused on the synthesis of ordered mesoporous alumina materials with specific surface area greater than 300 m2/g,ordered pore channels and narrow pore distribution.The ordered mesoporous alumina (OMA) material is expected to be widely used in the fields of catalysis,adsorption and separation(Yuan et al.,2008;Wu et al.,2011;Liu et al.,2016).At present,the research on application of ordered mesoporous alumina is still in the experimental stage and the synthesis method is still complex and difficult.Thus,it is of great significance to develop and study new ordered mesoporous alumina material.

It is generally believed that the interaction between the metal and the support plays an important role in sulfidation degree of the active component(Hu et al.,2020).The strong interaction between active metal Mo (W) and the traditional γ-Al2O3support would result in the low sulfidation degree of MoS2sulfide phase.Moreover,the metal Ni and γ-Al2O3is easy to form the non-active spinel phase,leading to the low utilization rate of metal Ni and the formation of undesired active type I "Ni-Mo-S" active phase (Van Veen et al.,1993).Therefore,optimizing the active-support interaction is the focus of researchers to improve the catalytic activity of the bimetallic catalysts in the HDS reaction.

Chelating agents,such as citric acid (CA) (Zhang et al.,2017;Li et al.,2011),ethylenediaminetetraacetic acid (EDTA) (Ortega-Domínguez et al.,2017;Badoga et al.,2012),aminotriacetic acid(NTA)(L′elias et al.,2009),and cyclohexanediaminetetraacetic acid(Cy-DTA) (Hiroshima et al.,1997) are commonly used as additives during the impregnation process to modify the active metal species.The hydrogenating catalysts prepared by using chelating agents have displayed excellent activity and stability,and the effects of complexing agents on the hydrogenating catalysts include the following three aspects:1) To delay the vulcanization of the auxiliary metal.2) To weaken the strong interaction between the active component and the support,thus contributing to form more highly active type II "Ni-Mo-S" active phases;3) To adjust the dispersions of the active phases on the surface of the catalyst support.The MoS2phase with a suitable sheet length and stacking number can effectively improve the HDS activity of the catalysts(Cattaneo et al.,2001;Nikulshin et al.,2014;Asadi et al.,2019).

Based on above information,this research aims at synthesizing ordered mesoporous alumina (OMA) with high specific surface area,concentrated pore size distribution and outstanding hydrodesulfurization performance.The NiMoE/OMA series catalysts were prepared by adding different amounts of chelating agent EDTA.Dibenzothiophene (DBT) was served as the model reactant to evaluate the hydrodesulfurization activity of NiMoE/OMA series catalysts.

2.Experimental section

2.1.Support preparation

Hydrothermal precipitation method is used to prepare ordered mesoporous alumina materials and the specific preparation steps are listed as follows:a certain amount of template PEG was added into the 0.6 mol/L aluminum nitrate solution.Then,1.2 mol/L ammonium carbonate solution was dropwise added to the above mixed solution at a stirring rate of 400 r/min.Continue to stir for 3 h and then age for 6 h in a 70°C water bath to obtain a white sol.After filtering,it needs to be dried at 70°C for 12 h and calcined in a N2atmosphere at 350°C for 3 h,then heated up to 550°C at a rate of 3°C/min for 6 h.Finally,the series ordered mesoporous alumina materials were obtained through changing the ratio of PEG/Al(0.05,0.1,0.2),which were named as OMA-0.05,OMA-0.1 and OMA-0.2,respectively.

2.2.Catalyst preparation

The supported NiMo catalyst precursor was prepared by incipient wetness impregnation method with the mass fraction loadings of 3.5 wt% NiO and 15 wt% MoO3respectively.NiMo/OMA catalyst precursor was obtained by dried at 80°C for 5 h and calcined at 550°C for 6 h.

The modified catalysts were prepared the same as above,except that different proportions of EDTA chelating agent are added during the impregnation process.The synthesized catalysts were named as NiMoE(X)/OMA,where X represents the ratio of EDTA/Ni.

2.3.Measurement and characterization

X-ray powder diffraction (XRD) was used to perform phase analysis on the samples by Shimadzu X-6000 X-ray powder diffractometer (Cu Kα radiation,40 kV tube pressure,and 30 mA tube current).Scanning range of 2θ is 0.5°-5°;wide angle range of 2θ is 5°-80°.

Fourier Infrared Spectroscopy (FTIR) was used to analyze the skeleton structure of the sample by the DIGILAB FTS-3000 Fourier Infrared Spectrometer of Tianmei Technology Company.The sample and KBr were mixed and compressed at a mass ratio of 1:100 with a resolution of 2 cm-1.Wavelength range:400-1200 cm-1.

Pyridine adsorption infrared spectrometer(Py-FTIR)was used to carry out qualitative and quantitative analysis of acidity in the samples by Digilab FT-IR from Bole Pacific Company.A vacuum(1 × 10-3Pa) was used to purify the sample for 2 h at 350°C,and then cooled to room temperature.

Scanning Electron Microscope (SEM) was used to analyze the micro-morphology of the samples by the Quanta 200F scanning electron microscope of the Dutch company FEI.

High resolution transmission electron microscopy(HRTEM)was used to calculate the dispersity,average length (Lav) and stack number (Nav) of the active phases on the sulfided catalysts by the JEM 2100 transmission electron microscope of Japan JEOL Company.The stacks number and average length of the MoS2active phases on each catalyst were statistically calculated by at least 300 sticks according to the equations (a) and (b) (L′opez-Benítez et al.,2017):

wherelirepresents the length of MoS2microstrips,niis the numbers of MoS2microstrips,Nirefers to the layers stack of MoS2active phases.

The fraction of Mo atoms located on the edges of MoS2clusters,which was denoted asfMo,was calculated using Mototal(the total number of Mo atoms)and Moedge(the number of Mo atoms located on the edges of MoS2particles).The values offMowere calculated by the equations (c) and (d).

WherefMorepresents the ratio of the number of Mo atoms at the edge to the total number of Mo atoms on the active phase,nirefers to the number of Mo atoms on the side of the MoS2crystal strips.

Low temperature N2-adsorption and desorption experiment(BET)was used to determine the specific surface area,pore volume and pore size distribution of the samples by using Micromeritics ASAP 2010 adsorption analyzer.The sample was pretreated at 100°C for 1 h under vacuum at a pressure of 15 μm Hg,and purified and degassed at 350°C for 3 h,and then subjected to static adsorption analysis at -196°C with N2.

Raman Spectrum was used to study the state of the metal species on the surface of the sample by using Renishaw's Raman spectrometer.The laser light source has a wavelength of 325 nm and a power of 8 MW.

X-ray Photoelectron Spectroscopy(XPS)was used to analyze the surface Mo species of the catalyst samples reduced by CS2.The XPSPEAK 4.1 software was used to Gaussian fit the Mo 3d XPS spectrum,and the relative content of Mo species in different valence states on the catalyst surface was calculated.The degree of sulfuration of Mo species and the utilization rate of Ni metal were calculated according to formula (e) and (f) respectively.

2.4.Catalytic activity measurements

1g of sieved 40-60 mesh catalyst is charged into the constant temperature section of the reactor.The two ends of the bed are filled with 20-40 mesh quartz sands,and each bed is separated by quartz cottons.3.0 wt% CS2cyclohexane solution was used to presulfide the oxidation state of the series catalysts.Firstly,the flow rate of the presulfiding solution was maintained at 60 mL h-1for 20 min,and then adjusted to 5 mL h-1.After 4 h of presulfiding,the device was washed three times with pure cyclohexane,and the remaining presulfiding solution was washed out.After pretreatment,a 500 ppm DBT cyclohexane solution was passed as the raw material under the conditions of the hydrogen pressure of 4 MPa,the temperature of 340°C,and the volume ratio of H2/Oil of 200.

3.Results

3.1.Characterization results of supports

3.1.1.XRD characterization

Fig.1(A) and (B) are the small-angle XRD and wide-angle XRD spectra of OMA materials prepared with different PEG/Al ratios respectively.As shown in Fig.1 (A),a relatively strong diffraction peak appears in the range of 0.5-1°,indicating that the pore structure of the synthesized alumina material is ordered.It can be seen from Fig.1(B) that there are no characteristic peaks of crystalline alumina,confirming that the synthesized alumina is an amorphous alumina material.

3.1.2.N2 adsorption-desorption characterization

Fig.1.Small-angle XRD a and wide-angle XRD b patterns of OMA materials with different PEG/Al ratios.

Fig.2 shows the pore size distribution (A) and N2adsorptiondesorption isotherm (B) of the alumina material under different PEG/Al ratios.It can be seen from Fig.2(A)that the pore sizes of all the synthesized alumina materials are concentrated in the range of 5-10 nm,proving the relatively ordered pore size distribution of the support.Fig.2(B) shows that the N2adsorption-desorption isotherm of the synthesized OMA material belongs to a characteristic type IV adsorption equilibrium curve with a H2 hysteresis loop in the range of P/P0=0.5～0.9.Table 1 summarizes the textural properties of OMA the series materials,among which the OMA-0.1 material shows an excellent pore structure with high specific surface area (328 m2g-1),large pore volume 0.74 (cm3·g-1) and big average pore size (8.1 nm).Therefore,OMA material with PEG/Al ratio of 0.1 is chosen as the ideal support for hydrodesulfurization catalysts.

3.2.Characterization results of catalysts

3.2.1.XRD characterization

Fig.3 shows the wide-angle XRD patterns of NiMoE/OMA catalysts modified with different EDTA/Ni molar ratios.As can be seen from Fig.3,NiMoE/OMA the series catalysts show a characteristic peak at 2θ=67.30°,indicating that the support exists in amorphous structure rather than crystalline alumina.Moreover,no peaks of metal oxides in the spectrum are detected,which indicates that there are no large particle aggregations of Mo and Ni species.

3.2.2.Raman characterization

Fig.4 shows the Raman diagrams of NiMoE/OMA oxidized state catalyst with different EDTA/Ni molar ratios.It can be seen from Fig.4 that the diffraction peaks at 339 cm-1and 850 cm-1are weak,which are attributed to the characteristic bending vibration peak of the Mo=O bond of MoO42-tetrahedral species.The MoO42-tetrahedral species are the products of the strong interaction between the support and the active metal,which are difficult to be vulcanized into the MoS2active phases in the pre-vulcanization stage(Xiao et al.,2018).In Fig.4,the series NiMoE/OMA catalysts show obvious wide peaks at 953 cm-1attributed to the stretching vibration of the Mo=O bond of two-dimensional polymer Mo7O264-species,which is easy to be vulcanized into MoS2active phase in the pre-vulcanization stage(Parola et al.,2002;Wang et al.,2015).Therefore,the Raman result shows that the NiMoE/OMA catalysts modified by EDTA possess an optimal support-metal interaction force.The characteristic peak of the modified NiMoE(1.0)/OMA catalyst at 953 cm-1is stronger than that of other catalysts,meaning that the MoS2active phases are easier to be formed on the NiMoE(1.0)/OMA catalyst when a moderate amount of EDTA is added.

3.2.3.Py-IR characterization

To analyze the acid types and acid amounts,the Py-TR spectra of NiMoE/OMA series catalyst is shown in Fig.S1.Fig.S1(A)is the Py-TR spectra at 200°C,representing the total acid amount,while Py-TR spectra at 350°C represent the strong and medium strong acids in Fig.S1(B)(Barzetti et al.,1996;Zhang et al.,2008).It can be found that the diffraction peak intensity of EDTA-modified NiMoE/OMA series catalysts increases firstly and then decreases followingly as the increase molar ratio of the EDTA/Ni.Table 2 shows the calculated acid amounts of NiMoE/OMA series catalysts.The order of the quantities of total acidities and the medium and strong acidities is NiMoE(1.0)/OMA ＞ NiMoE(1.5)/OMA ＞ NiMoE(0.5)/OMA ＞NiMoE(0.2)/OMA ＞NiMo/OMA.

3.2.4.XPS characterization

Fig.S2 shows the Mo 3d XPS spectra of the series sulfided NiMoE/OMA catalysts.According to the peak splitting data in Fig.S2,the binding energy peaks at 228.9 ± 0.1 eV and 232.0±0.1 eV are assigned to the Mo 3d5/2 and Mo 3d3/2 of Mo4+,and the binding energy peaks at 230.5±0.1 eV and 233.6±0.1 eV are ascribed to the Mo 3d5/2 and Mo 3d3/2 of Mo5+,whereas the binding energies of Mo 3d5/2 and Mo 3d3/2 of Mo6+are 232.5 ± 0.1 eV and 235.6 ± 0.1 eV,respectively (Jos′e et al.,2002).According to Table 3,the sulfidation degree of the NiMoE/OMA catalysts modified by EDTA is all higher than 0.6 compared to the unmodified NiMo/OMA catalyst.With the increasing molar ratios of EDTA/Ni,the sulfidation degrees of the corresponding catalysts are also increased (Li et al.,2019).

Fig.2.Pore diameter distribution a and N2 adsorption-desorption b patterns of OMA materials with different PEG/Al ratios.

Table 1 Textural properties of OMA materials with different PEG/Al ratios.

Fig.3.The XRD patterns of NiMoE/OMA oxided catalyst with different EDTA/Ni molar ratios.a NiMoE(0.2)/OMA,b NiMoE(0.5)/OMA,c NiMoE(1.0)/OMA,d NiMoE(1.5)/OMA.

Fig.S3 shows the spectra of Ni 2p XPS of the series sulfided catalysts.In Fig.S3,the existence forms of Ni metals in the support are mainly NiMoS,NiSxand NiO compounds,among which the corresponding binding energies are 856.2 ± 0.1eV,the 853.2±0.2eV and 861.2±0.2eV respectively(Lai et al.,2013).It can be found that Ni metals are mainly formed as NiMoS active phases,whereas the peak of NiSx is almost non-existent over NiMoE(0.5)/OMA and NiMoE(1.0)/OMA catalysts.As shown in Table 4,the proportion of NiMoS phase increases from 0.54 over NiMo/OMA to 0.72 over NiMoE(1.5)/OMA,indicating that the addition of EDTA could significantly promote the formation of NiMoS active phases.

Fig.4.Raman patterns of NiMoE/OMA oxided catalysts.a NiMo/OMA,b NiMoE(0.2)/OMA,c NiMoE(0.5)/OMA,d NiMoE(1.0)/OMA,e NiMoE(1.5)/OMA.

Table 2 The amounts of acid determined by Py-FTIR of series catalysts modified by EDTA.

3.2.5.HRTEM characterization

HRTEM images of NiMoE/OMA sulfided catalysts and the distributions of stacking layers of MoS2active phases are shown in Fig.5.The average length(Lav),average number of layers(Nav)and dispersion(fMo)of MoS2stacks are calculated and listed in Table 5.With the increase of EDTA/Ni ratios,the average lengths of NiMoE/OMA series catalysts decreases from 3.5 nm of NiMoE(0.2)/OMA to 3.2 nm of NiMoE(1.0)/OMA and then increases to 3.3 nm of NiMoE(1.5)/OMA.Among all the catalysts,the NiMoE(1.0)/OMA catalyst exhibits short average length (3.2 nm) and suitable stacking layers(2.8).ThefMovalue is closely related to the layer number and the average length of active phases,and the follows trend of NiMoE(1.0)/OMA＞NiMoE(0.5)/OMA＞NiMoE(0.2)/OMA ＞NiMoE(1.5)/OMA.

Table 3 The Mo 3d XPS analysis results of the sulfided NiMoE/OMA series catalysts.

Table 4 The Ni 2p XPS analysis results of the sulfided NiMoE/OMA catalysts.

Table 5 Average length (Lav) and average layer number (Nav) of MoS2 crystallites.

Fig.5.HRTEM images of the series sulfided NiMoE/OMA catalysts.a NiMoE(0.2)/OMA,b NiMoE(0.5)/OMA,c NiMoE(1.0)/OMA,d NiMoE(1.5)/OMA.

3.3.The results of DBT hydrodesulfurization

Fig.6 shows the DBT HDS performance of EDTA modified catalysts at different WHSVs.The HDS activity of NiMoE/OMA catalysts show a gradually increasing trend with the decrease of WHSVs from 100 h-1to 20 h-1.The DBT HDS conversions of the series catalysts are in order of NiMoE(1.0)/OMA ＞ NiMoE(1.5)/OMA ＞NiMoE(0.5)/OMA ＞NiMoE(0.2)/OMA,among which the highest desulfurization rate of NiMoE(1.0)/OMA catalyst is 97.7%at 20 h-1WHSV.Compared with the unmodified NiMo/OMA catalyst,NiMoE/OMA catalysts modified with EDTA display higher hydrodesulfurization activities at each WHSVs.

The cyclohexylbenzene (CHB) and biphenyl (BP) are the main products of DBT HDS detected from GC-MS,which are the main products of HYD and DDS reaction route respectively [24].The product distributions of the series NiMoE/OMA catalysts are shown in Fig.7.As can be seen,the CHB selectivity over NiMoE/OMA catalysts increases as the ratio of EDTA/Ni increases from 0.0 to 1.0,whereas the BP selectivity decreases as the increasing amount of EDTA.Therefore,the appropriate addition of EDTA has great influence on the promotion of HYD reaction route of DBT HDS to some extent.

4.Discussion

Through homogeneous precipitation method,the synthesized ordered mesoporous alumina (OMA) has ordered mesoporous channels(8.1 nm),the high surface area(328 m2g-1)and the large pore volume 0.74 (cm3·g-1),which could eliminate the diffusion resistance of the reactants and products through the mass transfer process and increase the accessibility of the reactants to the active metals.NiMoE/OMA catalysts synthesized by EDTA postmodification method show higher HDS activity of DBT compound compared with the unmodified NiMo/OMA catalyst.The HDS evaluation of DBT compound shows that as the molar ratios of EDTA/Ni increase,the DBT hydrodesulfurization conversions of the series catalysts increase at first and then decrease gradually,among which NiMoE(1.0)/OMA catalyst displays high DBT hydrodesulfurization activity,reaching 97.7%.The addition of EDTA has several following effects that are beneficial to the catalytic activity.

Fig.6.DBT HDS performance of NiMoE/OMA catalysts at different WHSVs.

Fig.7.Product distribution of the series NiMoE/OMA catalysts.

The physicochemical properties of the OMA support contributes a lot to the diffusion process of DBT molecules.Through adjusting the ratio of PEG/Al during preparing ordered mesoporous alumina materials,the OMA with PEG/Al ratio of 0.1 exhibits high surface area (328 m2g-1),large pore volume 0.74 (cm3·g-1) and big average pore size(8.1 nm).The well-ordered mesoporous structure and open pore channels of the OMA support can largely reduce the diffusion resistances of reactants and products compared to the traditional alumina materials.In addition,it can be seen from Fig.3,no peaks of metal oxides are detected in the XRD spectra,which indicates that the outstanding pore structure of OMA supports can also avoid the aggregation of active metals.

Secondly,the addition of EDTA is beneficial to the formation of Ni-Mo-S active phases.The XPS results show that NiMoE(1.0)/OMA and NiMoE(1.5)/OMA catalysts have the highest sulfidation degree of 0.65.The higher the sulfidation degree of the catalyst,the easier to form the type II “Ni-Mo-S′′active phase,which is more conducive to the HDS performance of the catalyst.As the molar ratio of EDTA/Ni increases,the proportion of NiMoS also increases,as shown in Table 4.It is well-known that NiMoS active phase is formed at the edge of MoS2layer structure by horizontal bonding(Liu et al.,2020).And some non-active phases are easily formed between Ni and alumina support,which leads to low utilization rate of Ni active metal.Coulier et al.,(2001)has reported that EDTA chelating ligand could stabilize nickel against sulfidation by forming a stable (NiEDTA)2-complex with the Ni promoter.The interaction between EDTA and Al cations acts as a driving force for decomposition of the(NiEDTA)2-complex so as to retard sulfidation of Ni to temperatures where Mo is completely sulfided in the form of MoS2(Al-Dalama,Stanislaus,2011).And Ni metallic species are combined with MoS2at the corner positions of hexagonal configurations to form the highly NiMoS active phase,which greatly enhance the utilization rate of Ni active metals.Therefore,the presence of EDTA chelating ligands can achieve a good level of interaction between Ni and the MoS2crystallites (Van Veen et al.,1993;Zepeda et al.,2016).

Moreover,the chelating agent EDTA promotes the dispersion of both Ni(II) and Mo(VI) on alumina.It has been reported that the chelating ligands EDTA tend to isolate Ni from the environment,thus avoiding the formation of excessive bulk Ni sulfide(Rana et al.,2007;Zhao et al.,2006).Meanwhile,EDTA had a higher coordinating constant with the surface Al cations than Mo(VI).Free EDTA can compete with Mo(VI) on the alumina surface sites.Therefore,chelating agent EDTA can also limits strong interaction between the metal Mo ions and alumina.According to the HRTEM images,The NiMoE/OMA-1.0 catalyst has a short length of MoS2active phase(3.2 nm),a moderate number of stacking layers (2.8),and a maximum dispersion degree (0.35).The Raman characterization results show that the stretching vibration of the Mo=O bond of polymerspecies (which is easy to vulcanize) over EDTAmodified NiMoE/OMA series catalysts is higher than that of unmodified NiMo/OMA.

The addition of EDTA can also change the acid properties of the modified catalysts so as to enhance DBT HDS activity.As shown in Py-FTIR characterization results,with the addition of EDTA increasing,the total acidities of NiMoE/OMA catalysts increase firstly and then decrease,among which the NiMoE(1.0)/OMA catalyst possesses the highest acid amount.It is well known that acidity plays an important role in improving the adsorption of DBT via its aromatic rings,thus finally enhancing HDS activity (Wang et al.,2016).The NiMoE(1.0)/OMA catalyst with the maximum total acid amount exhibits the high selectivity of CHB product,which indicates that the acid properties of the NiMoE/OMA catalysts have great influence on HYD and DDS pathways during the DBT HDS process.

Based on the above analysis,the introduction of an appropriate amount of EDTA can promote the HDS performance of DBT and HYD reaction route due to the combination effect of the weak interaction between the metal support,the high sulfidation degree of NiMoS phase and suitable acid properties on HDS.

5.Conclusions

Ordered mesoporous alumina (OMA) was successfully synthesized by homogeneous precipitation with aluminum nitrate as the inorganic aluminum source,ammonium carbonate as the precipitant and polyethylene glycol (PEG) as the template.The chelating agent EDTA was added to adjust the interaction between the support and the metal,so as to improve the sulfidation degree of Mo2S and enhance utilization rate of Ni metal.NiMoE/OMA catalysts synthesized by EDTA post-modification method displayed higher HDS activities of DBT compound compared with the unmodified NiMo/OMA.And among all the synthesized catalysts,NiMoE(1.0)/OMA catalyst showed the highest desulfurization rate of 97.7%and more HYD hydrogenation products of CHB.

Acknowledgments

This research is financially supported by the National Natural Science Foundation of China (No.21878330,21676298) and the National Key R&D Program of China (2019YFC1907602) and the CNPC Key Research Project(2016E-0707).

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.petsci.2021.11.005.