Diagnosis and Monitoring of Volatile Fatty Acids Production from Raw Cheese Whey by Multiscale Time-Series Analysis

: Anaerobic treatment is a viable alternative for the treatment of agro-industrial waste. Anaerobic digestion reduces organic load and produces volatile fatty acids (VFA), which are precursors of value-added products such as methane-rich biogas, biohydrogen, and biopolymers. Nowadays, there are no low-cost diagnosis and monitoring systems that analyze the dynamic behavior of key variables in real time, representing a signiﬁcant limitation for its practical implementation. In this work, the feasibility of using the multiscale analysis to diagnose and monitor the key variables in VFA production by anaerobic treatment of raw cheese whey is presented. First, experiments were carried out to evaluate the performance of the proposed methodology under different operating conditions. Then, experimental pH time series were analyzed using rescaled range ( R / S ) techniques. Time-series analysis shows that the anaerobic VFA production exhibits a multiscale behavior, identifying three characteristic regions (i.e., three values of Hurst exponent). In addition, the dynamic Hurst exponents show satisfactory correlations with the chemical oxygen demand (COD) consumption and VFA production. The multiscale analysis of pH time series is easy to implement and inexpensive. Hence, it could be used as a diagnosis and indirect monitoring system of key variables in the anaerobic treatment of raw cheese whey.


Introduction
Cheese whey is a waste generated in cheese production, which is characterized by its high organic load [1]. Different biological and physicochemical treatments have been proposed to recover compounds with added value, such as proteins and lactose. However, a significant number of small and medium industries do not have the necessary economic resources to implement technologies for the treatment of whey, which leads to the inadequate disposal of the waste, causing serious problems of environmental pollution [2]. Anaerobic treatment is an alternative for the valorization of this effluent. The anaerobic digestion of cheese whey reduces the organic load with simultaneous recovery of value-added by-products such as VFA, methane-rich biogas, biohydrogen, and biopolymers [3,4].
Anaerobic digestion (AD) is a complex and intrinsically unstable process. Indeed, the anaerobic treatment is formed by consecutive stages conducted by a specific group of microorganisms with different metabolic characteristics which are not fully understood. Moreover, variations of the input variables (hydraulic flow rate, influent organic load) and the effects of inhibitory agents may easily lead to the wash-out condition [5].
Thus, proper diagnosis and control of AD are of utmost relevance to avoid serious breakdowns and guarantee operational stability [6,7]. To this end, it is necessary to monitor the variables that provide information on the performance of each stage. The key variables of each stage can be defined by their relationship with the performance of the associated microbial consortium. Indeed, carbohydrate and protein degradation are related to hydrolysis, VFA production with acidogenesis-acetogenesis, and biogas production with acetogenesis-methanogenesis [8]. However, at the industrial level, in situ and in-line measurements in AD plants are limited to pH, temperature, water flow, biogas flow, level, and pressure [9][10][11].
Although pH and biogas production can be used as indicators of the overall status of the AD process, they fail to indicate the stress level of the reactor in real time [12]. On the one hand, biogas production can show the overall process performance. However, a decrease in biogas production often occurred after the process was severely inhibited. On the other hand, pH is not a suitable parameter for a well-buffered bioreactor since pH has low sensitivity. VFA concentration has been suggested as an ideal indicator [12]. Indeed, VFA accumulation is usually interpreted as the result of methanogenesis inhibition or organic overloading, leading to a risk of process breakdown. Thus, for AD process diagnosis, several authors have proposed different thresholds based on VFA measurements as early warning indicators to judge the reactor status [5,[13][14][15][16].
Hence, for the AD diagnosis, different model-based and data-driven estimation approaches have been proposed to determine key variables, including VFA concentration, using the available information in real time (e.g., temperature, pH, conductivity measurements). Model-based estimation methods are based on the reconstruction of unmeasured variables from the measurements available through the mathematical model of the process [11,17]. Model-based estimation approaches applied to AD processes include state-observers [18][19][20][21][22], Bayesian methods [23], and filtering techniques [24,25]. The most significant limitation of model-based approaches is that they depend on the accuracy of the process model, i.e., the lower the precision of the model, the greater the estimation error. Moreover, model uncertainties degrade model predictions, and the measurements are affected by noise and disturbances.
Data-driven methods are based on statistical models trained to learn trends from collected historical [26]. Data-driven estimation approaches are gaining broad interest at the industrial level because of their practicability, robustness, and flexibility to be developed and applied to a wide range of processes, in addition to their independence from a mathematical process model [17,27]. Examples of data-driven methods used in AD processes include artificial neural networks (ANN) [28,29], regression models [30], machine learning [16,31], and fuzzy inference systems [32]. Data-driven methods based on multiscale approaches have been introduced due to the superior robustness and accuracy concerning single-scale analysis [33]. Indeed, the motivation for using multiscale methods for feature extraction is that the data from almost all practical processes are of a multiscale nature. For instance, the events can occur at different locations and with different localization in time and frequency. Moreover, conventional statistical methods based on a single scale are not ideally suited for separating the deterministic component from the normal process variation. Lee et al. [34] applied a multiscale approach with principal component analysis and wavelets to a full-scale biological anaerobic process. Data-driven methods have two main disadvantages: (i) High implementation costs in terms of design effort (architecture selection, training stage, extensive multivariable data sets) and computational cost. (ii) The lack of phenomenological interpretability since they are based on black-box models to predict the desired estimated states.
In the last few years, multiscale data-driven methods based on fractal analysis have shown great potential to infer the key variables of the process. Méndez-Acosta et al. [35] proposed the rescaled range (R/S) analysis of pH time series obtained from an up-flow reactor for the anaerobic digestion of tequila vinasses, where high correlations between the multiscale indices and critical variables of the process (i.e., COD, VFA, and biogas production) were determined. Similar results were found by Hernandez-Martinez et al. [36] applying multiscale analysis in pH time series obtained from anaerobic treatment of tequila vinasses in continuous digesters (FBR and CSTR). They concluded that the R/S analysis could qualitatively estimate essential variables for the process evolution. Finally, Sanchez-Garcia et al. [37] studied the anaerobic digestion of cheese whey, determining that multiscale analysis can be used for monitoring the anaerobic treatment of complex substrates. Hence, the above works have shown that through multiscale analysis of time series, it is possible to dynamically characterize the key variables of the anaerobic digestion process using different substrates. Multiscale fractal analysis has two nice features for its practical applicability. The multiscale analysis allows a better understanding of the underlying process phenomena occurring at different scales. Furthermore, a reduction of the implementational costs is obtained (in terms of both design effort and computational cost) since it is based on using a single easy-to-measure and low-cost variables, such as pH and electrical conductivity.
In this paper, the application of the multiscale fractal analysis is presented for diagnosis purposes of the anaerobic treatment of raw cheese whey. It should be noted that previous contributions on the multiscale fractal analysis of anaerobic digestion have been applied to stable operating conditions favoring methane production. However, multiscale analysis in situations where the process is affected by inhibition or problems that do not favor methane production has not been evaluated. An efficient monitoring system must assess the dynamic behavior of the variables in a broad region of operation, especially in the region where the process tends to destabilize since, in these conditions, it is when action is required to control the process. Indeed, the need for a proper diagnosis and monitoring system is more significant when the process is in critical conditions. The multiscale fractal analysis is performed on pH time series of two different experimental conditions of pH and temperature. The time series are analyzed using R/S multiscale analysis. The results indicate that it is possible to track the key variables using dynamically calculated multiscale indices, which open the possibility of developing reliable diagnosis and monitoring systems of the anaerobic processes that contribute to the operational improvement of the process.

Materials and Methods
In this section, the experimental conditions used for the anaerobic treatment of raw cheese whey are presented. The operating conditions, as well as the offline and online measurements performed during the experimental runs, are also described.

Experimental Setup
The anaerobic fermentation of raw cheese whey was carried out in an anaerobic sequencing batch reactor (AnSBR) with an effective volume of 3.5 L. The sludge used as inoculum was obtained from a laboratory-scale anaerobic digester fed with raw cheese whey OLR = 3.6 g COD L −1 d −1 and operated at 33 • C, pH = 8.0 and HRT = 30 d. Experimental runs were performed using fresh raw cheese whey as a substrate collected from a community dairy located in Acajete Veracruz, México. The cheese whey was dried, milled, and screened in a 1 mm mesh. After that, it was stored at 4 • C. Table 1 shows the values obtained of the characterization of cheese whey and inoculum.
In order to evaluate the effect of initial pH and temperature on VFA production, four experimental essays to different conditions (pH = 5.5-T = 35 • C, pH = 5.5-T = 40 • C, pH = 7.5-T = 35 • C, and pH = 7.5-T = 40 • C) were carried. In experimental trials performed, the pH control was only established to avoid acidification of the system, i.e., the pH is not less than the initial pH. The pH control for the pH increase was not established, causing pH to increase, so the initial pH value was not maintained.

Analytical Methods
The raw cheese whey and inoculum were characterized by determining the concentration of fat, protein, lactose, total carbohydrates (TCH), VFA, total solids (TS), volatile suspended solids (VSS), as well as the COD and alkalinity. Standard techniques were used to determine COD [38] and VSS [39]. Total carbohydrates and protein concentration were determined according to Dubois et al. [40] and Bradford et al. [41], respectively. An Agilent gas chromatograph (model 7820 A) equipped with a capillary column and a flame ionization detector was used to determine the concentrations of short-chain organic acids. Conversion factors were used to pass VFA, TCH, VSS, and protein concentration to COD equivalent concentration. The biogas production was measured with a Gow Mac 580 series gas chromatograph equipped with an isothermal column oven. Helium and Nitrogen were used as carrier gas. A NI cRIO-9074 system from National Instrument and an M300 pH/ORP transmitter from Mettler-Toledo were used to acquire pH time series. The pH meter has an accuracy of ±0.02, capturing 1 datapoint per second in each experiment. The pH time series shows temporal variations, reflecting the dynamic changes of the interactions between the characteristic phenomena involved in the process ( Figure 1). For all experiments, fluctuations are observed around the initial pH, followed by an increase in pH to a maximum value, to reach an apparent steady state, where the pH stabilizes finally. It can be observed that the pH dynamics depend on the process operating conditions. In the experiments at initial pH = 5.5, the system showed an increase in pH at a time greater than 35 h, while for pH = 7.5, such increase occurred at 80 h.

Multiscale Time Series Analysis
In the literature, several methodologies are available for multiscale time-series analysis. However, due to its simplicity and easy implementation, rescaled range (R/S) [42] is widely applied for time-series analysis obtained from physical, chemical, and biological processes.
The R/S statistical measures the range of the deviations of the partial sums in a time series about its average, rescaled by the series standard deviation. For a time series where X Ns is the subsequence mean and σ S (N S ) is the sample standard deviation, which are defined as The R/S statistic follows a power law, (R/S) = aN H S where a is a constant and H is the Hurst exponent. A log-log plot of (R/S) as a function of N S ∈ (N S,min , N S,max ), gives a straight line with slope H. If the series data is independent (e.g., white-noise process), the plot is roughly a straight line with slope H = 0.5. If H > 0.5, the time series

Multiscale Time Series Analysis
In the literature, several methodologies are available for multiscale time-series analysis. However, due to its simplicity and easy implementation, rescaled range (R/S) [42] is widely applied for time-series analysis obtained from physical, chemical, and biological processes.
The R/S statistical measures the range of the deviations of the partial sums in a time series about its average, rescaled by the series standard deviation. For a time series ZN = (Zi) with length N, it is considered a subsequence XNS = (Xi) of length NS, where NS < N. The R/S statistic for XNS is calculated as where Ns X is the subsequence mean and

Multifractal Analysis
Multifractality is a valuable tool for explaining many patterns seen in nature. Particularly, multifractal analysis allows for the investigation of a mixture of fractal dimensions characterizing the inherent complexity in some data series. Based on the thoughts of Barabasi and Vicsek [45] and Katsuragi and Honjo [46], the multifractal analysis can be done through the calculation of the rescaled range through the q-norm of σs, that is, In this case, the R/S statistical is given by (R/S) q = (R/S)σ(N S )/σ q (N S ), and the average range is expected to follow the scaling behavior (R/S) q = aN 2H q S , where H q is the q-th Hurst exponent. If H q is constant for all q then the underlying phenomena are monofractal. A nontrivial dependence of H q vs. q indicates that the process is multifractal. It should be recalled that a multifractal system is a generalization of a fractal system in which a single exponent (the fractal dimension) is not enough to describe its dynamics. In general, multifractality also indicates the nonlinear nature of the mechanisms that generated the series [47,48]. To describe the degree of multifractality, we propose a multifractal index (I M ) defined as, where H(q) is the H value as a function of q.

Results and Discussion
In this section, the results obtained from applying the multiscale analysis on the pH time series collected in the VFA production process by the anaerobic treatment of raw cheese whey are presented. First, the effect of pH and temperature on the overall performance of the process (i.e., VFA production and COD consumption) is described. Subsequently, the correlations obtained between the multiscale analysis and the experimental data of the key variables of the process are shown. Finally, the application of multifractal analysis to identify dynamic changes in the process associated with the different stages of the process is shown.

Experimental Tests
Whey is considered a complex substrate that is constituted mainly by carbohydrates, proteins, and lipids. At the conditions evaluated in this work, it is observed that total carbohydrates are consumed up to 95% in the first 35 h of operation ( Figure 2). However, the COD consumption only reaches a consumption between 30 and 50%, which suggests that the consumption of lipids and proteins was less significant. At pH = 7.5 and T = 40 °C, 50% of the COD consumption is reached, while at pH = 5.5 and T = 35 °C the minimum consumption is obtained with 30% ( Figure 3). These results correspond to those obtained by Perna et al. [49], who achieved COD removal percentages between 14% and 32% using diluted whey. However, the COD consumption only reaches a consumption between 30 and 50%, which suggests that the consumption of lipids and proteins was less significant. At pH = 7.5 and T = 40 • C, 50% of the COD consumption is reached, while at pH = 5.5 and T = 35 • C the minimum consumption is obtained with 30% ( Figure 3). These results correspond to those obtained by Perna et al. [49], who achieved COD removal percentages between 14% and 32% using diluted whey.
Regarding VFA production, the amount and types of VFA produced mainly depend on the initial pH. In Figure 3, the production of VFA is shown, noting that at initial pH = 5.5, butyric, valeric, and caproic acids are favored. In contrast, at initial pH = 7.5, the production of acetic and propionic acids is favored. These results agree with that reported by different authors, who indicate that at pH < 6, butyric acids are favored, and at pH > 6 propionic and acetic acids are favored [50][51][52][53].
The pH time series exhibit dynamic changes that correspond to the observed changes in key process variables. For the first 35 h, where the carbohydrate consumption rate is high, the pH series shows strong fluctuations, which could indicate the high activity of hydrolytic microorganisms. Subsequently, after 40 h, the experiments at pH 5.5 present an alkalinization stage, observing an increase in pH values of 6.8-7. This behavior can be associated with the production of biogas. Indeed, the presence of CO 2 promotes the generation of bicarbonates which causes an increase in pH. For the experiments at pH = 7.5, the alkalinization stage begins at 80 h, showing a slight increase in pH, reaching a maximum value of pH = 8.3. Finally, when the pH stabilizes, the highest VFA production rate is observed, which indicates the high activity of the acidogenic-acetogenic microorganisms. However, the COD consumption only reaches a consumption between 30 and 50%, which suggests that the consumption of lipids and proteins was less significant. At pH = 7.5 and T = 40 °C, 50% of the COD consumption is reached, while at pH = 5.5 and T = 35 °C the minimum consumption is obtained with 30% ( Figure 3). These results correspond to those obtained by Perna et al. [49], who achieved COD removal percentages between 14% and 32% using diluted whey.

R/S Analysis
R/S analysis is applied to the pH series of the four experimental conditions. The time scale interval is given by N S,min = 20 data (20 s) and N S,max = 3600 data (1 h), which according to Sanchez-Garcia et al. [37], is the time scale where the characteristic phenomena corresponding to each of the stages of the anaerobic digestion of crude whey can be identified. Figure 4 shows the statistic (R/S)2 as a function of the time scale Ns, observing that the power-law ((R/S) = aNs H ) exhibits three changes in slope (i.e., three different values of the Hurst exponent, H 1 , H 2 , and H 3 ), named as Zones 1, 2, and 3. The presence of different values of the Hurst exponent suggests that the R/S analysis can identify three characteristic phenomena immersed in the fluctuations of the pH series. These zones correspond to those reported in the anaerobic digestion of whey [37] and tequila vinasse [35,36], where each region was attributed to a different anaerobic digestion stage. Based on the above observations, Zone 1 is related to hydrolysis, Zone 2 corresponds with the acidogenesis and acetogenesis stages, and Zone 3 is related to methanogenesis. It should be noted that in conventional anaerobic digestion processes, the experimental conditions favor methane-rich biogas production. In contrast, in this study, the operating conditions favor the VFA production, so the methanogenic stage is not significant. each region was attributed to a different anaerobic digestion stage. Based on the above observations, Zone 1 is related to hydrolysis, Zone 2 corresponds with the acidogenesis and acetogenesis stages, and Zone 3 is related to methanogenesis. It should be noted that in conventional anaerobic digestion processes, the experimental conditions favor methane-rich biogas production. In contrast, in this study, the operating conditions favor the VFA production, so the methanogenic stage is not significant.  Table 2 shows the values of the Hurst exponent for each zone, where it is observed that the H 1 values are around 1.0, indicating that in this time scale, there are characteristic phenomena with similar behavior for all the operating conditions evaluated. On the other hand, the Hurst exponent values in Zones 2 and 3 show more significant variability between operating conditions, which indicates that the phenomena on these characteristic scales are more susceptible to changes in pH and temperature. According to the above, Zone 1 could be correlated with the hydrolytic stage. In all experiments, the COD consumption is mainly to carbohydrate degradation, which is consumed up to 95% in the first two days. Zones 2 and 3 could be correlated with the acidogenic-acetogenic stages, which show the formation of different fatty acids depending on the operating conditions evaluated. For example, at pH = 5.5, higher production of butyric and caproic acids is obtained, while at pH = 7.5, higher production of acetic and propionic acids is obtained. Note that the H 3 value is higher when there is a more significant formation of butyric and caproic acids.
In order to identify correlations between the key variables and the Hurst exponents, the R/S analysis is implemented considering 6 h windows, which are shifted every hour, covering the whole series. This allows one to obtain the Hurst exponent values as a function of scale and time. For the four experimental conditions, Figure 5 shows the dynamic R/S analysis, where the three zones described in Figure 4 can be identified, noting that the H values show temporal variations in each zone, which can be correlated with those of the process stages.
In order to avoid different interpretations of the dynamic Hurst exponent, the average is carried out in each zone, which leads to an H dynamic profile for each zone (i.e., Hd 1 , Hd 2, and Hd 3 ). These profiles are compared with the experimental measurements of COD and VFA.  Considering the four experimental conditions, Figure 6 shows the comparison between the COD measurements with Hd 1 , where it is observed that the Hd 1 dynamic exhibits the same trend observed with the COD consumption, which confirms that the time scale of Zone 1 corresponds to the hydrolysis stage in the anaerobic treatment of whey.   Figure 7 shows the comparison between the total VFA measurements produced with the Hd 3 dynamic profile, finding that both profiles exhibit the same behavior, suggesting that the time scale of Zone 3 corresponds to the acidogenic-acetogenic stage. Hd 2 profile shows behavior that cannot be correlated with any of the process variables.
Appl. Sci. 2021, 11, x 11 of 1 Figure 7 shows the comparison between the total VFA measurements produced wit the Hd3 dynamic profile, finding that both profiles exhibit the same behavior, suggestin that the time scale of Zone 3 corresponds to the acidogenic-acetogenic stage. Hd2 profil shows behavior that cannot be correlated with any of the process variables. Moreover, due to the operating conditions evaluated, there is no significant methan production, so the methanogenic stage is not considered as an active stage of the proces Then, Hurst exponents dynamic could be used as qualitative indexes for monitoring th COD and VFA for different operating conditions.

Multifractal Analysis
The multifractal analysis is a useful tool to identify information at different tim scales associated with the inherent complexity in bioprocesses.
The anaerobic treatmen of raw cheese whey is a highly complex process involving transport phenomena and bio reactions, affecting overall process stability. This interaction depends on the evaluate operating conditions and is reflected in the type and quantity of products generated.
In Figure 8, the application of the multifractal R/S analysis is presented, considerin q ∈ [0. 2,20]. For the experimental test at T = 40 °C and initial pH = 5.5, Figure 8a shows th statistic (R/S)q, noting that the slopes (Hq) of the three zones show variations at differen values of q-norm. Figure 8b shows that the Hq values of the three zones present nonlinea variations with respect to the q-norm that exhibit significant differences in each zone, sug gesting the presence of a multifractal behavior in the global process, as well as in each o the stages that make up the process. Moreover, due to the operating conditions evaluated, there is no significant methane production, so the methanogenic stage is not considered as an active stage of the process. Then, Hurst exponents dynamic could be used as qualitative indexes for monitoring the COD and VFA for different operating conditions.

Multifractal Analysis
The multifractal analysis is a useful tool to identify information at different time scales associated with the inherent complexity in bioprocesses.
The anaerobic treatment of raw cheese whey is a highly complex process involving transport phenomena and bioreactions, affecting overall process stability. This interaction depends on the evaluated operating conditions and is reflected in the type and quantity of products generated.
In Figure 8, the application of the multifractal R/S analysis is presented, considering q ∈ [0. 2,20]. For the experimental test at T = 40 • C and initial pH = 5.5, Figure 8a shows the statistic (R/S)q, noting that the slopes (H q ) of the three zones show variations at different values of q-norm. Figure 8b shows that the H q values of the three zones present nonlinear variations with respect to the q-norm that exhibit significant differences in each zone, suggesting the presence of a multifractal behavior in the global process, as well as in each of the stages that make up the process. The degree of multifractality of each process stage is evaluated with the multifracta index (IM). The higher the IM value, the greater the multifractal behavior. Table 3 shows the multifractal index of the three zones for the four experimental tests. The lowest values of the multifractal index correspond to Zone 1, associated with the hydrolytic stage. Re gardless of the operating conditions evaluated, a similar behavior is observed in all exper iments. Whereas Zone 3, associated with the acidic stage, shows the highest values of the multifractal index, which correspond to the generation of products obtained in this stage It should be noted that different types of VFA are produced and the quantity produced depends on the operating conditions. Globally, the multifractal index identifies which stage of the process presents more significant interactions between the underlying phenomena. In order to determine the changes in multifractality as the process is carried out, the multifractal index is calculated considering 6 h moving windows. Figure 9 shows the scale-temporal changes of the mul tifractal index, where it is observed that the dynamic changes of the IM correspond to the variations observed in the consumption/production rates of the key variables. In Zone 1 it is observed that the most important changes in IM occur in the first 30 h of the process where the maximum consumption of carbohydrates is carried out. In Zone 3, more signif icant dynamic changes are observed associated with the changes in the VFA production rate. Increases in the IM correspond to the increases in the VFA production rate. The degree of multifractality of each process stage is evaluated with the multifractal index (IM). The higher the IM value, the greater the multifractal behavior. Table 3 shows the multifractal index of the three zones for the four experimental tests. The lowest values of the multifractal index correspond to Zone 1, associated with the hydrolytic stage. Regardless of the operating conditions evaluated, a similar behavior is observed in all experiments. Whereas Zone 3, associated with the acidic stage, shows the highest values of the multifractal index, which correspond to the generation of products obtained in this stage. It should be noted that different types of VFA are produced and the quantity produced depends on the operating conditions. Globally, the multifractal index identifies which stage of the process presents more significant interactions between the underlying phenomena. In order to determine the changes in multifractality as the process is carried out, the multifractal index is calculated considering 6 h moving windows. Figure 9 shows the scale-temporal changes of the multifractal index, where it is observed that the dynamic changes of the IM correspond to the variations observed in the consumption/production rates of the key variables. In Zone 1, it is observed that the most important changes in IM occur in the first 30 h of the process, where the maximum consumption of carbohydrates is carried out. In Zone 3, more significant dynamic changes are observed associated with the changes in the VFA production rate. Increases in the IM correspond to the increases in the VFA production rate. Appl. Sci. 2021, 11, x 13 of 16

Conclusions
This work proposed the application of fractal analysis to pH time series collected from the anaerobic treatment of raw cheese whey in an anaerobic batch digester for four different operating conditions variating the pH and temperature. The rescaled range analysis identifies three regions that suggest three underlying phenomena in the anaerobic digestion treatment for raw cheese whey. The dynamic R/S analysis has permitted a qualitative description of two key variables for the four experimental sets. Thus, the R/S fractal methodology can be used as a tool for indirect monitoring of COD consumption and VFA production to different operating conditions. The fact that key variables can be monitored by pH measurements under different operating conditions is a promising advance for developing an economic diagnosis and monitoring and control system that can indicate process performance. Indeed, results in this paper can be used in conjunction with available diagnosis methods using thresholds based on VFA trends. Besides, the proposed multifractal index provides information on dynamic changes in consumption/production rates of key process variables, which are directly associated with the inherent complexity of the process.

Conclusions
This work proposed the application of fractal analysis to pH time series collected from the anaerobic treatment of raw cheese whey in an anaerobic batch digester for four different operating conditions variating the pH and temperature. The rescaled range analysis identifies three regions that suggest three underlying phenomena in the anaerobic digestion treatment for raw cheese whey. The dynamic R/S analysis has permitted a qualitative description of two key variables for the four experimental sets. Thus, the R/S fractal methodology can be used as a tool for indirect monitoring of COD consumption and VFA production to different operating conditions. The fact that key variables can be monitored by pH measurements under different operating conditions is a promising advance for developing an economic diagnosis and monitoring and control system that can indicate process performance. Indeed, results in this paper can be used in conjunction with available diagnosis methods using thresholds based on VFA trends. Besides, the proposed multifractal index provides information on dynamic changes in consumption/production rates of key process variables, which are directly associated with the inherent complexity of the process.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could appear to have influenced the work reported in this paper.