Identification of a mathematical model and parameter estimation of erythromycin migration in two different porous media , based on column tests

The amount of pharmaceuticals found in groundwater has risen over the last few years. Erythromycin is an example of an antibiotic widely used in human health care and veterinary practice that can be transported into the subsurface. The aim of the presented research was to 1) determine the mathematical model of erythromycin migration in two different porous media, 2) estimate the model parameters and 3) compare the migration of the antibiotic in the investigated media. The research was conducted in a specially prepared laboratory stand where column tests were performed. One column was filled with glass granules (70% SiO2, 600–800 μm in diameter) and the second column was filled with a natural sediment (sandur sand). The migration of a conservative tracer and erythromycin was examined in both columns. The experiments were performed in two separate steps. In the first step, a conservative tracer, subject to advection and dispersion processes, was injected into the coluand its transport was investigated. The second step involved investigating the migration of erythromycin. A conductivity meter was installed at the output of the column in order to determine tracer concentrations based on calibration curves. Short-time pulse injections and continuous injection were applied during the experiments. An interpretation of the experiments results was conducted in the MATLAB environment. A mathematical model of erythromycin migration was determined from the shape of pulse breakthrough curves, which were characterized by a set of descriptors: the time of maximum tracer concentration at the output tmax, the spread of the breakthrough curve s, and relative tracer recovery e. This procedure involved implementing an identification algorithm developed by the authors. It was proved that the migration of erythromycin is best described by a hybrid model that assumes the coexistence of equilibrium and non-equilibrium sorption. In the next stage of the research, the transport and sorption parameters were estimated through numerical optimization procedures. The convergence between theoretical and experimental breakthrough curves was analysed qualitatively by calculating the root mean square error RMSE and correlation coefficient r. Pharmaceuticals entered into the circulation of substances con­ nected with the hydrological cycle (Fig. 1) (KÜMMERER, 2008; MOMPELAT et al., 2009; ZAJĄC et al., 2013). The list of the pharmaceuticals which need to be monitored includes erythro­ mycin – a macrolide antibiotic used in human health care and veterinary practice (DE VOOGT et al., 2009; MA et al., 2015). Article history: Manuscript received June 16, 2017 Revised manuscript accepted January 24, 2018 Available online June 21, 2018


INTRODUCTION
Advances in medicine and the common use of pharmaceuticals have resulted in the presence of such compounds in both surface and groundwater as well as in geological sediments (WATKIN-SON et al., 2009;ZUCCATO et al., 2010;BU et al., 2013;WU et al., 2014;BOROŃ & PAWLAS, 2015;YAO et al., 2017). Pharma- Monika Okońska*

Abstract
The amount of pharmaceuticals found in groundwater has risen over the last few years. Erythromycin is an example of an antibiotic widely used in human health care and veterinary practice that can be transported into the subsurface. The aim of the presented research was to 1) determine the mathematical model of erythromycin migration in two different porous media, 2) estimate the model parameters and 3) compare the migration of the antibiotic in the investigated media. The research was conducted in a specially prepared laboratory stand where column tests were performed. One column was filled with glass granules (70% SiO 2 , 600-800 µm in diameter) and the second column was filled with a natural sediment (sandur sand). The migration of a conservative tracer and erythromycin was examined in both columns. The experiments were performed in two separate steps. In the first step, a conservative tracer, subject to advection and dispersion processes, was injected into the column and its transport was investigated. The second step involved investigating the migration of erythromycin. A conductivity meter was installed at the output of the column in order to determine tracer concentrations based on calibration curves. Short-time pulse injections and continuous injection were applied during the experiments. An interpretation of the experiments results was conducted in the MATLAB environment. A mathematical model of erythromycin migration was determined from the shape of pulse breakthrough curves, which were characterized by a set of descriptors: the time of maximum tracer concentration at the output t max , the spread of the breakthrough curve s, and relative tracer recovery e. This procedure involved implementing an identification algorithm developed by the authors. It was proved that the migration of erythromycin is best described by a hybrid model that assumes the coexistence of equilibrium and non-equilibrium sorption. In the next stage of the research, the transport and sorption parameters were estimated through numerical optimization procedures. The convergence between theoretical and experimental breakthrough curves was analysed qualitatively by calculating the root mean square error RMSE and correlation coefficient r.
ceuticals entered into the circulation of substances connected with the hydrological cycle ( Fig. 1) (KÜMMERER ed., 2008;MOMPELAT et al., 2009;ZAJĄC et al., 2013). The list of the pharmaceuticals which need to be monitored includes erythromycin -a macrolide antibiotic used in human health care and veterinary practice (DE VOOGT et al., 2009;MA et al., 2015). A R T I C L E I N P R E S S According to research findings, the effectiveness of removing erythromycin from the contaminated water in the sewage treatment process varies, and the percentage elimination of erythromycin depends on the technologies applied, i.e. the conventional or advanced treatment (WATKINSON et al., 2007;BOLEDA et al., 2011;CZECH, 2012;AL QARNI et al., 2016). Hence, there is a possibility that the antibiotic reaches the subsurface environment. According to the literature, cases of erythromycin presence in groundwater have already been reported (PENG et al., 2014;MA et al., 2015;YAO et al., 2017). Therefore, it is very important to investigate the conditions of erythromycin migration in natural porous media. BU et al. (2013), among others, pointed out the need for better understanding of the transport and fate of pharmaceuticals, including erythromycin in various media. Such study was conducted to recognize the leaching of erythromycin through porous media (SIEMENS et al., 2010), erythromycin adsorption (SUN et al., 2009) and the kinetic interaction with bacterial ribosomes (LOVMAR et al., 2004). The need for and the advantages of column experiments in the research of organic compounds transport were stressed by BANZHAF & HEBIG (2016).
The aim of this pilot study was to: 1) determine the mathematical model of erythromycin migration by means of the developed algorithm of model selection; 2) estimate the model parameters and 3) compare conditions of erythromycin migration in two different porous media. The authors used the MATLAB environment which enables a numerical solution of migration equations and an estimation of the parameters of the multiparameter models to be determined by means of numerical optimization procedures. The authors expect that the findings of the research can be further applied to mathematical modelling of erythromycin migration processes in the groundwater environment, including the construction of prognostic models.

MATERIALS AND METHODS
The test of erythromycin migration was conducted in a laboratory stand equipped with a filtration column 50 cm in length and 9 cm in diameter. Column tests were conducted for two different porous media. In the first test, the column was filled with spherical glass granules (70% SiO 2 , 600-800 µm in diameter). The porous medium was uniformly grained, with the coefficient of uniformity being U=1.20. The coefficient U is defined as: where D 10 is the grain diameter [mm] corresponding to 10% of the sample passing by weight (which means that 10% of the par-ticles are smaller than the diameter D 10 ) and D 60 is the grain diameter [mm] at 60% passing (HOLTZ & KOVACS, 1981).
The second experiment included a test using natural sediment, i.e. a sample of sandur sand from the area of Sierosław in the central Wielkopolska region, Poland. The investigated sand (reference Si-W1/1) was classified as medium quartz sand, with a grain uniformity coefficient of U=4.13 (OKOŃSKA, 2016). This value means that the sandur sand is less uniform in terms of grain size, compared to the glass granules. The total organic carbon content in the sand determined by the elemental analysis method was 0.114%. The sand sample contained trace amounts of clay. The clay minerals determined in the soil mineral composition by means of X-ray analysis were vermiculite, kaolinite and illite.
Distilled water was used to saturate the column with porous material and to prepare solutions containing tracers. During the column test, a sodium chloride solution was injected initially. This substance was regarded as a conservative tracer (CTR) (LEIBUNDGUT et al., 2009), i.e. an ideal tracer that undergoes the advection-dispersion processes. Next, a solution containing erythromycin (ERY), the migration of which was being investigated, was injected. It was a pulse tracer injection of known concentration ( Fig. 2A). In the next step -a continuous tracer injection was performed in accordance with Heaviside's step function ( Fig. 2B), (OKOŃSKA, 2006). Before the erythromycin injection, the conservative tracer was fully removed from the porous material. The concentration of the substance in the injected solution was C 0 =500 mg·dm -3 for sodium chloride and C 0 =1541 mg·dm -3 for erythromycin. The presence of erythromycin in the aquatic environment has been widely reported at lower concentrations. However, due to the methodology of the research and the size of the filtration columns the authors decided to impact the investigated medium with a strong identification signal so as to obtain a noticeable reaction of the system to that signal (SÖDERSTÖM & STOICA, 1988).
Electrolytic conductivity of water was investigated at the output of the filtration column, using a portable meter HQ40D with a laboratory graphite conductivity electrode INTELLICAL. The values recorded were then converted into tracer concentrations using the previously devised calibration curves.
During the tests, the total porosity coefficient n [-] was estimated by means of the weight method and the effective porosity coefficient n e [-] by means of the volumetric method. Water flow rate Q [cm 3 ·s -1 ] through a soil sample was also recorded. It allowed determination of the hydraulic conductivity k [m·s -1 ].
The mathematical models describing tracer migration through porous materials neglected chemical and bio-chemical reactions, biodegradation, chemo-and electro-osmotic and capillary effects. Moreover, it was assumed that the porous material cannot be deformed and is fully saturated with liquid solution and all the material and geometrical parameters present in the model are constant.
As for the conservative tracer, its transport through the material was described with a one-dimensional equation that takes account of the advection-dispersion processes (A-D model) (WE-BER et al., 1991): where C is the tracer concentration in the liquid phase [mg·dm -3 ], k is the hydraulic conductivity [m·s -1 ], i is the hydraulic gradient In the case of erythromycin, the migration mathematical model was identified with impulse breakthrough curve shape indicators (the so-called descriptors) and the authors' own model selection algorithm. The analysis included the following descriptors: time of maximum tracer concentration at the output t max [s]: where C max is the maximum tracer concentration at the output [mg·dm -3 ], spread of the breakthrough curve s [s]: relative tracer recovery e [%]: where t in is time interval of the tracer injection, when C(t) = C 0 at the input [s].
The model selection algorithm takes into account the percentage deviations of the descriptors that describe the erythromycin breakthrough curve with respect to descriptors determined for the conservative tracer breakthrough curve: D ERY is the descriptor value for the breakthrough curve of erythromycin, and D CTR is the descriptor value for the breakthrough curve of a conservative tracer.
Seven of the eleven considered sorption models were selected for analysis: • four simple models, taking account of equilibrium sorption or non-equilibrium sorption, are described by the following equation (WEBER et al., 1991): where the variable s was substituted with an equation describing one of the following sorption models: the Henry model (H model), the Freundlich model (F model), the Langmuir model (L model), or the reversible kinetic sorption model (R model); three hybrid models that assume the co-occurrence of equilibrium sorption (s e ) and non-equilibrium sorption (s n ), described by a system of equations (CAMERON & KLUTE, 1977): where variables s e i s n are substituted with equations describing one of the following models: the Henry model with reversible sorption (H-R model) the Freundlich model with reversible sorption (F-R model) the Langmuir model with reversible sorption (L-R model) where r s is the density of the porous medium [kg·dm -3 ], n is total porosity [-], K H is the Henry distribution coefficient [dm 3 ·kg -1 ], k 2 is the first reversible sorption rate coefficient [dm 3 ·kg -1 ·s -1 ], k 3 is the second reversible sorption rate coefficient [s -1 ], K F is the Freundlich sorption coefficient [dm 3 ·kg -1 ], n F is the Freundlich sorption exponent [-], a L is the Langmuir constant [dm 3 ·kg -1 ], and b L is the total sorption capacity of the solid phase [mg·kg -1 ]. The next stage of the research involved a value identification (estimation) of transport and sorption parameters in the MAT-LAB environment based on the recorded breakthrough curves. The estimation procedure included: • numerical solutions of partial differential equations of transport and sorption assuming suitable boundary conditions, • the application of the numerical optimization method covering both the global and local optimization procedure (MARCINIAK et al., 2009).
The equations were solved by means of the pdepe solver which utilizes the numerical finite element method (FEM). For identification calculations, the authors used the Optimization Toolbox and the in-built lsqcurvefit function for the non-linear solution of the optimization problem. This function, featuring the large-scale optimization option, utilizes the area-restricted search method (COLEMAN & LI, 1994).
The model was entered into with initial values of the investigated parameters and with column test data. In the advectiondispersion model the parameters investigated were the value of hydraulic coefficient k and the value of longitudinal dispersivity a. The subject of identification in sorption models was the sorption parameters. In those models, the value of the hydraulic coefficient k and the value of longitudinal dispersivity a were substituted with parameter values determined from the conservative tracer tests.

A R T I C L E I N P R E S S
Searching the point where the error function assumes the smallest value in the search range went on until the set optimization criteria were met. While assessing the quality of the fitting of the theoretical tracer breakthrough curve to the experimental breakthrough curve the root mean square error RMSE and correlation coefficient r were used: where A is the set of estimated parameters, C m is the concentration value [mg·dm -3 ] measured in time t m [s], C t is the concentration value obtained in the simulation [mg·dm -3 ]; where cov signifies the covariance between the values observed C m [mg·dm -3 ] and calculated C t [mg·dm -3 ], and SD is standard deviation.
The applied procedure of parameter estimation in the MAT-LAB environment was verified by determining the dispersion parameters in the FIELD programme that uses the analytical solutions of the advection-dispersion model and the optimization method based on the least squares method (MALOSZEWSKI, 1981;HENDRY et al., 1999;FOHRMANN et al., 2001). In this programme the fitting of the theoretical curve to experimental data is determined by a model efficiency indicator ME [%], where ME = r·100.

RESULTS
The breakthrough curves of the conservative tracer (CTR) and erythromycin (ERY) migration through glass granules are shown in Fig. 3A, and the breakthrough curves of the conservative tracer and erythromycin migration through sand are shown in Fig. 3B. In the second case, the continuous curve for the conservative tracer has not been registered due to the anticipated long time of the column test, which could cause clogging of the deposit. Figure 3 shows erythromycin migration retardation compared to the conservative tracer. In the diagrams, time t of substance migration through the porous media refers to time t 0 that corresponds to distance x [m] and average water flow velocity v [m·s -1 ]. Measured concentration C is presented relative to the injected concentration C 0 .
The value of the hydraulic coefficient k, that was determined in the column test (Table 1) and entered into the advection-dispersion model as the initial value, in the identification procedure in the MATLAB environment changed only within the set measurement error (± 10-15%). As a result of the interpretation of the experimental breakthrough curve of the conservative tracer, the following values of longitudinal dispersivity were obtained: a=0.0007 m for the glass medium and a= 0.0096 m for sand (Table 2). Figure 4. shows the theoretical conservative tracer breakthrough curves calculated in the MATLAB environment against the background of experimental points. A very good degree of breakthrough curve fitting was obtained: the RMSE ranging between 0.014-0.036 and the correlation coefficient r=1 ( Table 2). The calculations made in the FIELD programme for impulse curves confirmed the determined values of the longitudinal dispersivity a (0.0009 m and 0.0088 m for glass granules and sand respectively, for model efficiency 99.8% and 98.9%).

A R T I C L E I N P R E S S
Results obtained for the two different porous media showed that longitudinal dispersivity differed by an order of magnitude, despite the same distance covered by the migration tracer. This would imply the dependence of parameter a on the degree of heterogeneity of the porous media (Fig. 5A).
For the advection-dispersion model it was assumed that the share of molecular diffusion in the conservative tracer migration through the media is negligible, for glass granules in particular. The parametric analysis of the impact of molecular diffusion coefficient D M [m 2 ·s -1 ] on the result of the calculations of the advec-    (Fig. 5B). Mass recovery of erythromycin was 99.1% (± 5%) for the glass medium, and 96.4% (± 5%) for the sand. On the basis of calculated deviations of descriptors, according to the identification algorithm, three mathematical models were selected for erythromycin and the analysed glass medium for model tests: R, H-R and F-R, and two models were identified for erythromycin and sand: L and L-R.
The estimation of migration parameters was carried out for the five indicated sorption models plus two simple models: H and F (Table 3). In the case of both porous media the best results were obtained for hybrid models which combine equilibrium and nonequilibrium sorption. In three cases, the best fit of theoretical curves to experimental data was obtained for the Henry model with reversible sorption (H-R model). In the case of sand and the continuous curve it was the Freundlich model with reversible sorption (F-R model). In the process of optimization the exponent n F =1 was often obtained for the Freundlich model. Therefore, this equation often became a linear isotherm equation. Figure 6. and Figure 7. show selected theoretical erythromycin breakthrough curves calculated in the MATLAB environment. A very accurate fit of the theoretical curves with experimental data was obtained. In the case of the H-R model, the RMSE falls within the range of 0.002-0.008, with a correlation coefficient r=0.992-0.999 (Table 3). It can be assumed that all of the considered hybrid models characterise erythromycin migration through the investigated porous media types very well.
The value of retardation R, estimated from the relationship t 0ERY to t 0CTR for impulse breakthrough curves, is 1.1 for glass medium and 1.4 for sand. It implies a low level of sorption intensity of erythromycin (WITCZAK et al., 2013).

CONCLUSIONS
The research conducted into erythromycin migration through the two types of porous material resulted in the identification of a mathematical model for describing antibiotic migration and esti-mation of the parameters present in the migration model. It can be stated that erythromycin migration is characterised accurately by hybrid models that assume both equilibrium sorption and nonequilibrium sorption. Best results were obtained for the Henry model with reversible sorption, while for the continuous breakthrough curve of erythromycin migration through sand -for the Freundlich model with reversible sorption.
The descriptors and the identification procedure correctly identified mathematical models for the glass medium. In the case of natural soils and a higher value of uniformity coefficient, only the type of sorption -reversible or irreversible -can be indicated. The identification algorithm was developed on the basis of sensitivity studies conducted by one of the authors in the past (MAR-CINIAK et al., 2013). These studies described the relationship between the curve's shape and the parameters of mathematical models for column tests performed on a uniformly grained porous material.
The calculated values of Henry distribution coefficients K H from the impulse curve are: 1.93E-02 dm 3 ·kg -1 and 1.77E-02 dm 3 ·kg -1 for the glass medium and 7.41E-02 dm 3 ·kg -1 for the sand. The value of the Freundlich sorption coefficient K F determined for sand and erythromycin from the continuous curve is 2.00E-02 dm 3 ·kg -1 . The Freundlich sorption exponent n F equal 1.88 may indicate, according APPELO & POSTMA (1999), "unfavorable exchange" which means that the fraction in the solution is larger than the fraction adsorbed. The first reversible sorption rate coefficient k 2 and the second reversible sorption rate coefficient k 3 are within the range of 1.76E-05 to 2.24E-06 dm 3 ·kg -1 ·s -1 and 2.72E-03 to 4.86E-06 s -1 for both porous media. The obtained retardation R values imply a low level of intensity of erythromycin sorption in the case of the investigated porous media (slightly higher for sand).
Making a comparison of the obtained sorption parameter values with the results of research conducted by other authors is difficult due to the variation in conditions under which such research was performed. The predicted distribution coefficient log-K OC for erythromycin is 2.53 dm 3 ·kg -1 (from American Chemical Society after SIEMENS et al., 2010). The authors received similar results for sand (logK OC =1.81 dm 3 ·kg -1 ), converting the Henry    , 1999). For glass granules the coefficient estimated from the impulse curve was: logK OC =0.29 dm 3 ·kg -1 . The obtained results (i.e. the identified mathematical models) confirm the findings made by BANZAF & HEBIG (2016) who claim that column tests are suitable for determining solute transport as they allow the determination of sorption under dynamic conditions, and they can provide the resulting non-equilibrium reactions.
The observed relationships between dispersivity a and the coefficient of uniformity U are consistent with the study of XU & ECKSTEIN (1997). The dispersivity value increases with the heterogeneity of porous media.