Kaili Wang

Guanhua Huang

Impact of Permeability on Flow and Solute Transport in Highly Heterogeneous Aquifer

Keywords: permeability, transport, Monte-Carlo, heterogeneous, moments

The impact of permeability on solute transport process in a highly heterogeneous aquifer was analyzed using the Monte-Carlo method. The logarithm of the permeability (lnK) of aquifer was considered as a non-stationary field with increments being a truncated fractional Lévy motion (fLm) generated using the SRA3DC code. MODFLOW and MT3DMS code were used to solve the flow and solute equations, respectively. Results indicate that larger C lead to a more heterogeneous permeability field. As C increases, solute plumes show more significant anomaly with sharper leading edge and wider tailing edge. The solute plume extends and its second spatial moments increase as C increases, while the first spatial moments of the solute plume are independent of C values. The longitudinal macrodispersivity is scale-dependent and increases as a power law function of time. Increasing C results in an increase in longitudinal macrodispersivity. It was concluded that the nature of transport process is highly dependent on the heterogeneity of the permeability field.

Introduction

Groundwater contamination has become one of the most important environmental issues all over the world. It is necessary to predict flow and contaminant spreading in the subsurface for the control of groundwater quality. However, the heterogeneity of porous media and incomplete knowledge of data information lead to difficulty in the estimation of hydraulic properties and geophysical variables, and thus, bring about difficulty in estimating or predicting subsurface flow and transport. The stochastic methods have been developed and used to deal with these difficulties [START_REF] Dagan | Flow and Transport in Porous Formations[END_REF]. There are mainly two categories of stochastic methods in widespread use including the moment equation method (MEM) and the Monte-Carlo method (MCM). Compared with the MEM, the MCM is conceptually and computationally straightforward, and most importantly, once the properties' distributions are established in details, the MCM can be applied to highly heterogeneous media. Up to present, the MCM has been widely used in researches on flow and transport process with the more accurate subsurface simulation models and the availability of high-speed supercomputers [START_REF] Hassan | Monte Carlo studies of flow and transport in fractal conductivity fields: Comparison with stochastic perturbation theory[END_REF][START_REF] Huang | Monte Carlo study of conservative transport in heterogeneous dual-porosity media[END_REF][START_REF] Molz | A fractal-based stochastic interpolation scheme in subsurface hydrology[END_REF][START_REF] Yan | Impacts of the spatial variation of permeability on the transport fate of solute plume[END_REF][START_REF] Painter | Prediction uncertainty for tracer migration in random heterogeneities with multifractal character[END_REF].

In the MCM simulation, the first step is to determine the distributions of the properties of geological formations. Researches have shown that the permeability (K) has a relative large spatial variability with its maximum value several orders being larger than its minimum value, while the spatial variability of porosity and other geophysical parameters are relatively small [START_REF] Neuman | A quasi-linear theory of Non-Fickian and Fickian subsurface dispersion.1. Theoretical analysis with application to isotropic media[END_REF]. In classical stochastic theories, K is treated as a lognormal, stationary and spatially correlated random function, which can be characterized by its mean, variance and single finite-range integral scale [START_REF] Dagan | Flow and Transport in Porous Formations[END_REF]. It has been proved that such a treatment can be only used for those media with relatively small heterogeneities [START_REF] Painter | Random fractal model s of heterogeneity: the Lévy-stable approach[END_REF]. In recent years, several researches have shown that the geological media can be highly heterogeneous and exhibit very large and abrupt changes in K field [START_REF] Liu | Multifractal analyses of hydraulic conductivity distributions[END_REF][START_REF] Painter | Evidence for non-Gaussian scaling behavior in heterogeneous sedimentary formations[END_REF]. And such a hydraulic field with continuously evolving heterogeneities can be described using the fractal scaling model [4, 11, and 12]. One of the commonly used fractal model to characterize the non-stationary K field is the fractional Brownian motion (fBm), with which the probability density function (PDF) of lnK increments for any given lag are assumed to be Gaussian. However, researches show that the lnK increments of heterogeneous formations are often non-Gaussian with an increased peak value around its mean and heavy (slow-decaying) tails [START_REF] Painter | Evidence for non-Gaussian scaling behavior in heterogeneous sedimentary formations[END_REF].

These non-Gaussian behaviors can be characterized using the probability distributions of the Lévy-stable family, which is referred to as the fractional Lévy motion (fLm).

Compared with the Gaussian distribution, the Lévy distribution can better mimic the sharp property contrasts associated to the geological formations. However, the heavy tails also give the Lévy-stable distribution an infinite variance. Painter [START_REF] Painter | Stochastic interpolation of aquifer properties using fractional Lévy motion[END_REF] argued that the increments cannot have an infinite-variance distribution because the incremental values must be bounded in order to prevent porosity or permeability from taking on nonphysical values, and strictly speaking, this is correct. The truncated PDF have finite variance with power-law decay tails over a large but finite range, resulting in behavior similar to that of the infinite-variance PDF. Therefore, the truncated PDF for lnK increments will be used in the following analysis.

Several approaches have been developed to generate fBm or truncated fLm fields.

These methods include the midpoint displacements method, the successive random addition (SRA) method [START_REF] Peitgen | The science of fractal images[END_REF], the Fourier filtering method [START_REF] Peitgen | The science of fractal images[END_REF], and the modified turning band method [START_REF] Yin | New methods for simulation of fractional Brownian motion[END_REF]. Lu et al. [START_REF] Lu | An efficient, three-dimensional, anisotropic, fractional Brownian motion and truncated fractional Lévy motion simulation algorithm based on successive random additions[END_REF] recently compared these methods in detail and concluded that with the SRA algorithm it is easy to understand the geometry and the scaling properties of the stochastic fractals and it is associated with a very attractive computational speed. They developed a simple and efficient three-dimensional simulation code named SRA3DC in Fortran to generate the fBm and the truncated fLm K fields. In addition, dispersional and quantile analyses were employed to analyze the increments of the generated fBm and truncated fLm fields to validate the algorithm.

Many researches have been conducted to investigate solute transport in porous media with fBm and / or fLm permeability fields. For example, Painter [START_REF] Painter | Prediction uncertainty for tracer migration in random heterogeneities with multifractal character[END_REF] investigated the statistical moments of travel time for a conservative tracer in the media with an fBm K field and a bounded fLm K field. Herrick et al. [START_REF] Herrick | Hydraulic conductivity, velocity, and the order of the fractional dispersion derivative in a highly heterogeneous system[END_REF] investigated the tritium transport at the highly heterogeneous MADE site; they found that the heavy-tailed K fields gave rise to the heavy-tailed velocity fields. Kohlbecker et al. [START_REF] Kohlbecker | Heavy-tailed log permeability distributions imply heavy-tailed log velocity distributions[END_REF] found that a heavy-tailed Levy-stable lnK increment field resulted in a heavy-tailed increment field of logarithm velocity (ln v). However, to our knowledge, few researches were about the analysis of plume spatial moment and dispersivity in those media with fLm permeability fields.

In this study, we will use the MCM to investigate solute transport process in the highly heterogeneous aquifer with a truncated fLm field for the lnK increments. The SRA3DC code will be modified and used to generate lnK fields for the MCM analysis. MODFLOW2000 code [START_REF] Mcdonald | A modular three-dimensional finite-difference ground water flow model[END_REF] and MT3DMS code [START_REF] Zheng | MT3DMS: A modular three-dimensional multi-species transport model for simulation of advection, dispersion and chemical reactions of contaminants in ground water systems: Documentation and user's guide[END_REF] will be used to solve the flow and solute transport problems, respectively. The simulated concentration will be averaged to obtain the ensemble mean concentration for analyzing the plume spatial moments and the longitudinal macro-dispersivity.

Theoretical Consideration

Lévy-stable Distributions

The Lévy-stable distributions have been studied since the 1930s [START_REF] Lévy | Théorie de l ' Addition des Variables Aléatoires[END_REF]. Their probability density functions (PDFs) do not have a closed analytic form except for some special cases. Instead, they are expressed by their characteristic function:

} )} 2 tan( ) ( 1 { exp{ ) ( iuk k sign i k C k + - - = απ β φ α α 1 ≠ α . ( 1
)
where α is the Lévy index with a range of (0,2], β is the skewness parameter with a range of [-1,1], C is the width parameter in the interval (0,∞), μ is the location parameter, k is the Fourier variable, i 2 =-1, tan() represents the tangent function and sign(k) is the sign function of the variable k. These characterizing parameters α, C, β and μ describe the index of stability, the spread, the skewness, and the location of the density, respectively .

When the Lévy-stable distribution is a standard symmetric function with β and μ values equal to zero, Eq. ( 1) can be simplified as:

α α φ k C k - = exp( ) ( . (2)
The PDF of symmetric Lévy-stable distribution can be calculated by inverse Fourier transform of Eq. ( 2)

∫ ∞ - = 0 )} cos( exp{ 1 ) ( dk kx Ck x L α α π . (3)
In Eq. ( 3) the Lévy index α distinguishes the PDFs from one another and characterizes the degree of deviation from the Gaussian distribution which corresponds to the special case

) ( 2 x L = α
. The width parameter C characterizes the discrete extent of variables relative to the mean value, which is similar to the standard deviation of a Gaussian distribution.

Fractional Lévy Motion

The fractional Lévy motion (fLm) is defined as a random and continuous function which has Lévy-stable increments over any lag. The statistical self-affinity or scaling property can be expressed as [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes[END_REF]

h H rh H h rh C r C r C C = ⇒ = α α α . ( 4
)
where C rh and C h are the width parameters for lag distance rh and h, respectively. H is the Hurst coefficient, quantifying the degree of interdependence between the increments of fractal process.

When α=2, the Lévy-stable distribution becomes a Gaussian distribution, and the fLm model is the same as the fBm model. rh σ and h σ are the standard deviation of the Gaussian distribution for lag rh and h, respectively, corresponding to the width parameters of the Lévy-stable distribution. Substituting rh σ and h σ into Eq. ( 4) yields:

h H rh H h rh r r σ σ σ σ = ⇒ = 2 2 2
.

(

) 5
In general, the Lévy-stable distribution of increments and the scaling property determine the essential features of fLm process. In this paper, we modified Lu's SRA3DC code based on these two features to generate a two-dimensional lnK field characterized by the truncated fLm distribution. It is necessary to point out that the width parameters C 1 and C 2 are replaced in the modified code by

α α 1 ) 2 1 ( ) 2 1 ( 2 2 2 0 1 - - = H H C C . ( 6
) 2 1 2 ) 2 1 ( H C C = . ( 7
)
where C 0 is the given width parameter for random numbers at corner points, C 1 and C 2 are the width parameters for random numbers at center points and edge middle points, respectively (see Lu et al. [START_REF] Lu | An efficient, three-dimensional, anisotropic, fractional Brownian motion and truncated fractional Lévy motion simulation algorithm based on successive random additions[END_REF] for more information).

A truncated fLm realization with H=0.3, C=0.2 and α=1.3 was generated for lnK field, with which quantile analysis was used to analyze the increments of lnK to estimate C and α. The relationship between the estimated C and the lag is shown in

Numerical Simulation

Problem Setup

A two-dimensional confined aquifer with a domain size 128 m x 128 m was considered with the grid size equaling to 1m (as shown in Fig. 2). Upstream and downstream boundaries were specified as constant head with a head difference of 0.7m in the mean flow direction. No-flow conditions were prescribed at the transverse, upper and lower boundaries. The porosity was considered to be spatially homogeneous with a value of unity for convenience. The geometric mean of lnK was chosen to be 2.5m/d and the increments of lnK were treated as a truncated fLm distribution. The local longitudinal and transverse dispersivities were 1m and 0.1m, respectively. To avoid the boundary effects, a total mass of 1kg was released at (10m, 64m) as a small point source of contamination. And the duration time was finished on the first transport step size.

Monte-Carlo Simulation

The MCM was used to simulate solute transport process in the domain. As shown in Table 1, three width parameters, i.e. C=0.05 (case 1), 0.1(case 2) and 0.2 (case 3) were chosen to analyze the impact of width parameters on solute transport process, respectively. Five hundred realizations for each case were carried out in the MCM to satisfy the convergence requirements. The simulation algorithm consists of the following steps:

(1) Using the modified SRA3DC code to generate the random lnK field with its increments following the truncated fLm distribution of given random seed number, mean, C, H and α values, respectively. The corresponding K can be obtained with the exponential transformation of the generated lnK random field for each realization.

(2) Solving the flow equation to obtain the velocity field for each realization of the K field. Firstly, flow problem was solved using MODFLOW with the Pre Conditioned Gradient solver (PCG2). Secondly, the cell by cell flow terms were then extracted from the MODFLOW solution and converted from mass flux to velocity.

(3) Solving the transport equation to obtain the concentration field for each realization of the velocity field using MT3DMS with the central finite difference scheme for advection solution. Each realization was stopped when 10% of the initial mass exited the domain.

(4) Repeating steps (1)-( 3) for all the realizations and averaging over all realizations to obtain the ensemble mean concentration distribution. The average concentration distribution will be used for moment analysis

Moment Analysis

Spatial moments of solute concentration are widely used to analyze solute transport process. The zero-order spatial moment indicating the mass of the system can be expressed as [START_REF] Zhang | Applied geostatistics in environmental sciences[END_REF]:

∫∫ Ω = dxdy t y x c t M ) , , ( ) ( 0 θ . (8)
where M 0 (t) is the mass of the solute in the system at time t; θ is the effective porosity; c(x, y, t) is the concentration at location (x, y) for time t, and Ω represents the interest area.

The first-order moments about the origin are the mass center of solute plume, which can be expressed as:

∫∫ Ω = xdxdy t y x c M t x c ) , , ( 1 ) ( 0
θ . ( 9a
) ∫∫ Ω = ydxdy t y x c M t y c ) , , ( 1 ) ( 0 θ . (9b)
where x c (t), y c (t) are the first moments representing the centroid coordinate of solute plume along the longitudinal and transverse directions at time t, respectively.

The second-order moments about the plume's center of mass characterizing the macrodispersion can be expressed as:

2 2 0 2 0 ) ( ) , , ( 1 )
)( , , ( 1 )
( t x dxdy x t y x c M dxdy x x t y x c M t M c c xx - = - = ∫∫ ∫∫ Ω Ω θ θ . ( 10
) 2 2 0 ) ( ) , , ( 1 )
- = ∫∫ Ω θ . ( 11
)
where M xx (t), M yy (t) are the second central moments of the concentration plume along the longitudinal and transverse directions at time t, respectively.

The longitudinal macrodispersion coefficient D can be expressed as follows:

dt dM D xx 2 1 = . ( 12
)
Under the assumption of a uniform flow field, the longitudinal macrodispersivity is defined as

dt v dM t xx x 2 1 ) ( = α . ( 13
)
where v is the average seepage velocity along the x coordinate, which can be determined by Darcy's law

θ h K v ∇ = . ( 14
)
where h ∇ is the hydraulic gradient, K is the geometric mean of the permeability, θ is the porosity which can be treated as a constant.

Results and Discussions

Permeability Field Properties

Velocity Field Properties

Fig. 4 shows the contours of velocity magnitude (v x 2 +v y 2 ) 1/2 in the domain for a realization of case 1-3. It can be found that the increase in C leads to larger variation in the velocity field. The largest value of the velocity magnitude is about 0.12, 0.24 and 0.45 for case 1, 2 and 3, respectively. These results are consistent with the permeability field properties. This is attributed to Darcy's law on condition that the hydraulic gradient and porosity are constant in this study. Consequently, the differences in the permeability field leads to the same trend of the differences in the velocity field. Combined with the research of Kohlbecker et al. [START_REF] Kohlbecker | Heavy-tailed log permeability distributions imply heavy-tailed log velocity distributions[END_REF] and our present study, it was necessary to make further investigation of the probability distribution of the velocity field. In this study, we mainly examine the impact of lnK with the increments following the truncated fLm distribution on the resulting increments of lnv. For convenience, the longitudinal component velocity vector (i.e. v x ) is chosen to be discussed. Fig. 5 is the sampled probability density plot of the increments for lag=1 in lnv x for for different C (i.e. case 2 and 3) values. It indicates that the probability densities of the increments in lnv x deviate significantly from Gaussian distribution and shift from the body to the tails, which can be approximated well by the Levy-stable distribution. These results are consistent with the results of Kohlbecker et al. that a heavy-tailed distribution of increments in lnK results in a heavy-tailed distribution of increments in lnv. Furthermore, as C increases, the probability densities cover a greater range of values, signifying a lager variation of v x . This is consistent with the variations in the velocity magnitude as show in Figure 4. In next section, we will discuss solute transport process obtained from these heavy-tailed velocity fields.

Contour Maps of Solute Concentration

Contour maps of solute concentration are one of the best ways to interpret solute transport process. Fig. 6 shows the comparison of contour maps for case 1-3 at 400 days including a homogeneous case (i.e. C=0) for better analysis. It can be found that solute plume is in Gaussian distribution for C=0 in Fig. 6(a). As C increases, it cover a grater range with a corresponding decrease in peak value, i.e., 7-8 for case 1 in Fig. As can be seen from Eq. ( 14), the mean velocity is a linear function of the mean of K with constant hydraulic gradient and porosity. Hence, it was concluded that the first spatial moments are determined by the mean of lnK and independent of the width parameter. Fig. 8 shows the second spatial longitudinal and transverse moments for different C (i.e. case 2 and 3) values. Fig. 9 shows the longitudinal macrodispersivity α x obtained using Eq. ( 13) with the second longitudinal moments from Fig. 8. It's known that for the hypothetical case of homogeneous conditions, with constant velocity and dispersion coefficients and instantaneous point injection, solute plume is Gaussian and the second moments increase linearly with time while the longitudinal macrodispersivity is constant Our results show that M xx and M yy increase linearly with time on the double logarithmic scale, implying that the second spatial moments can be approximated by a power law function of time, which leads to α x increase as a power law function of time as shown in Fig. 9. These results may be attributable to that the hydraulic properties are random variables characterized by the Lévy-stable distribution, resulting in the spatial distribution of solute concentration deviates from Gaussian distribution and the second spatial moments of solute plume are nonlinear functions of time. This implies that non-Ficken dispersion occurs. Zou et al. [START_REF] Zou | Analytical solutions to non-Fickian subsurface dispersion in uniform groundwater flow[END_REF] have discussed the parabolic relationship between the dispersivity and time, and they have argued that the non-linearly dispersive process is caused by the heterogeneity of the hydraulic properties. Additionally, M xx , M yy and α x increase as C increases, meaning that larger degree of heterogeneity causes larger dispersion of solute plume.

Conclusion

In this study, we used the Monte-Carlo simulation method to investigate the sensitivity of solute transport to the characteristics of the highly heterogeneous and non-stationary random fields. Based on the comparison results as discussed in Section 4, the following conclusions are reached:

Firstly, the heterogeneity of the permeability field is sensitive to the values of the fLm characterizing parameters. A more heterogeneous permeability field is characterized by a larger C value. This is determined by the features of C, which characterizes the discrete extent of variables relative to the mean value and a larger C value leads to a larger variation in K fields.

Secondly, a larger C value value lead to a larger variation in the velocity field which are consistent with the permeability field properties. The investigation of the probability distribution shows that the probability density plot of the increments in lnv x can be approximated well by the Lévy-stable distribution. In addition, as C increases, the probability densities cover a greater range of values, signifying a larger variation of velocity.

Thirdly, contour maps of solute concentration show that larger C leads to more significantly anomalous transport with a sharper leading edge and a wider tailing edge in the plume. The basic reason for this phenomenon is that larger C lead to heavy-tailed K distribution, which directly give rise to heavy-tailed velocity field and then results in more significant anomalous transport.

Fourthly, the first spatial moments of the solute plumes are independent of the fLm characterizing parameters C. The second spatial moments can be approximated by a power law function of time, which leads to the longitudinal macrodispersivities increase as a power law function of time. In addition, both the second spatial moments and the longitudinal macrodispersivities increase as C increases, meaning that a larger degree of heterogeneity causes a larger dispersion of solute plume.

The results summarized above imply that solute transport process is highly dependent on the heterogeneity of the permeability field. The sensitivity analysis underscores the need for careful aquifer characterization for accurate estimation or prediction of subsurface flow and transport. However, our present study only considered two-dimensional case with relatively simple boundary conditions. Further investigation will be conducted for three-dimensional case in highly heterogeneous aquifer with more complex conditions.

Fig. 1 (Fig. 1 (

 11 Fig.1(a). The fitted H value is 0.3259 which is close to the input 0.3 value, and the

Fig. 2 .

 2 Fig. 2. Schematic illustration of the computational domain

Fig. 3

 3 Fig. 3 represents realizations of the lnK fields for case 1-3. It can be found that the spatial distribution patterns of the lnK fields are similar for different C values.However, the increase in C values leads to the wider disperse range of the lnK field.Specially, the value of lnK deviates from 2.5 gradually represented by the enlarged

Fig. 3 .

 3 Fig. 3. Comparison of the lnK field for (1) case 1; (2) case 2; (3) case 3

Fig. 4 .

 4 Fig. 4. Comparison of the contours of velocity magnitude for (1) case 1; (2) case 2; (3) case 3.

Fig. 5 .

 5 Fig. 5. Comparison of probability density plot of the increments for lag=1 in lnv x for case 2 and 3

  6(b), 5-6 for case 2 in Fig. 6(c) and 3-4 for case 3 in Fig. 6(d), respectively. Compared with the homogeneous case, solute plumes for larger C values have a more significant anomalous shape with a sharper leading edge and a wider tailing edge. The concentration gradient decreases with the increase of the C value. All these features reflect more significantly anomalous transport process in a more highly heterogeneous media.

Fig. 6 .

 6 Fig. 6. The contour maps of the solute concentration for different width parameters C at 400 days: (a) C=0 (i.e. homogeneous case); (b) case 1; (c) case 2; (d) case 3

Fig. 7 .

 7 Fig. 7. The first spatial moments (i.e. x c and y c ) of the solute plume versus time for case 2 and 3

Fig. 8 .Fig. 9 .

 89 Fig. 8. The second spatial moments of the solute plume versus time for case 2 and 3: (a) Moments at the longitudinal direction M xx ; (b) Moments at the transverse direction M yy

Table 1 .

 1 Chatacerizing parameyters considered in this study

	Case	α	C 0	H
	1	1.3	0.05	0.3
	2	1.3	0.1	0.3
	3	1.3	0.2	0.3


Huishan Lu

Zhengguang An

Huanyu Jiang
email: hyjiang@zju.edu.cn

Yibin Ying
email: ybying@zju.edu.cn

Discrimination between mature vinegars of different geographical origins by NIRS

The feasibility of near infrared spectroscopy (NIRS) for discrimination between mature vinegar of different geographical origins (Taiyuan and Qingxu, China) is presented in this research. NIR spectra were collected in transmission mode in the wavelength range of 800-2500 nm. Qualitative analysis models were developed based on Principal Component Analysis (PCA) and Discriminant Analysis (DA). The prediction performance of calibration models in different wavelength range was also investigated. The best models gave a 96.3% classification of mature vinegars of the two geographical origins in the range of 800-2500 nm. The results demonstrated that NIRS technique could be used as a rapid method for classification of geographical origin of mature vinegars.

Introduction

With increasing concerns about food products adulteration, falsification of added-value claims and traditional foods protection, food authentication is an essential challenge that must be faced in many different quality control tasks, such as guaranteeing the actual origin of a product and detecting deliberate or accidental adulteration of valuable food components by inferior ingredients. The food industry, regulatory authorities and consumers are all interested in authenticating raw materials and verifying the origin or the specific character of a product in order to promote quality and to assure food safety.

In particular, assurance of the quality and origin of vinegars has become an issue of supreme commercial importance as a means for searching for safety in the vinegar industry and for fighting to curb unfair competition, also taking into account the huge diversity of vinegars available nowadays on the market and the differences in the final sale price depending on a wide range of quality factors. Vinegar is a product resulting from the alcoholic and subsequent acetic fermentation of any fermentable starting material rich in carbohydrates. The final quality and organoleptic properties of vinegar are determined by a profound chemical complexity conditioned by both the raw material employed and the particular elaboration process used. Traditional and selected vinegars are produced by slow acetification methods, which usually require long periods of maturation in wood to reach the desired acetic degree. These traditional elaboration processes, involving simultaneous acetification and ageing, produce excellent quality vinegar, but at the expense of increasing production time and cost, which is used to justify their premium price. Therefore, reliable methods are required for quality and economical reasons in order to objectively differentiate vinegars according to their origin and commercial type, and thus provide real protection to superior quality products and to ensure vinegar authenticity. Flavour and aroma are some of the most important factors in the determination of vinegar character and quality. A large variety of compounds covering a wide range of volatilities and concentrations are responsible for the sensorial complexity of vinegars. The volatile fraction of any specific vinegar, and thus its aromatic composition, is significantly influenced by many factors, the main sources of variation being the particular raw material used, the acetification system employed for its production, and eventually wood ageing. Therefore, since the volatile profile of vinegar represents a fingerprint of the sample, it is reasonable to assume that a classification approach based on the analysis of volatile components would be an efficient tool for evaluating vinegar authenticity. In fact, a number of previous studies have described the effective application of the analysis of volatile compounds for the characterization and differentiation of different vinegars. Although it should be noted that this is only a feasibility study, the promising results obtained justify a similar approach to be considered in future in order to better evaluate its actual performance and to broaden the field of application to a wider range of vinegar types.

Materials and method

Samples

In this research, a total of two geographical origins of mature vinegars were obtained in different region. Thirty six bottles of mature vinegar samples of "Donghu" brand were from Qingxu, and thirty six bottles of "Ninghuafu" brand from Taiyuan. In all, 12 samples were of 1 year age, 12 samples were of 3 year age, 12 samples were of 5 year age in "Donghu" Group and "Ninghuafu" Group, respectively. Before the experiment, the mature vinegar samples were stored in the laboratory at a constant temperature of 25±1 0 C for more than 48h to have an equalization room temperature. The samples were all original vinegar liquid without dilution. Sixty mature vinegar samples were used in the calibration set, whereas, 12 samples (2 for each age of each group) were selected as the validation set.

Spectral measurements

Samples taken from freshly opened bottles of mature vinegar were scanned in transmission mode using a commercial spectrometer Nexus FT-NIR (Thermo Nicolet Corporation, Madison, WI,USA) which was equipped with an interferometer, an InGaAs detector, and a broad band light source (Quartz Tungsten Halogen, 50 W). Samples were scanned in a 1 mm optical path-length rectangular quartz cuvette with air as reference at room temperature. NIR spectra were collected using OMNIC software (Thermo Nicolet Corporation, Madison, WI, USA) and stored in absorbance format. The spectral range was from 800 to 2500 nm, the mirror velocity was 0.9494 cm s-1, and the resolution was 16 cm-1 in this work. The spectrum of each sample was the average of 32 successive scans.

Chemometrics and data analysis

Chemometrics analyses including PCA and DA were performed using a commercial software package (TQ Analyst software, Thermo Electron Corp., Madison, Wisc.). Spectra were exported from OMNIC software in absorbance format to TQ Analyst software before being analyzed.

Principal Component Analysis (PCA)

In this work, PCA was used to reduce the dimensionality of the NIR transmission spectra recorded on mature vinegar samples and extract the most useful and relevant information. PCA transforms the original independent variables (wavelengths) into new axes and calculates the principal components (known as PCs) as new variables by relevant algorithm to replace the original data [START_REF] Mouazen | Classification of soils into different moisture content levels based on NIS-NIR spectra[END_REF]. The PCs are orthogonal, so the data presented on the new axes are uncorrelated with each other [START_REF] Xu | Methods of Chemometrics[END_REF]). The PCs account as much as possible for the variability in the original variables [START_REF] Cozzolino | Feasibility study on the use of visible and near-infrared spectroscopy together with chemometrics to discriminate between commercial white wines of different varietal origins[END_REF]. Each spectrum will have its own unique set of scores; therefore, a spectrum can be represented by its PCA scores in the factor space instead of by intensity in the wavelength space [START_REF] Park | Near-infrared diffuse reflectance for quantitative and qualitative measurement of soluble solids and firmness of delicious and gala apples[END_REF]. The scores of such PCs replace the original spectra data to be used for further modeling.

Discriminant Analysis (DA)

DA was used in this work to build a classification model to discriminate "Ninghuafu", and "Donghu" mature vinegar samples based on predefined classes and determine the percentage of correct classifications. Because the first ten PCs cover most of the variation (>99.9% of the total variance) of the raw spectral data, the first ten PCs attained through PCA calculation of the spectral data were employed for DA. In TQ Analyst, discrimination of the groups consists of calculating the Mahalanobis distance of a sample from the centers of gravity of the considered groups; one can then clarify the properties that distinguish the different groups. A sample is assigned to the group that has the shortest Mahalanobis distance to it, with which it is considered to be the most similar.

Results and discussion

spectral analysis

Figure 1 shows the average NIR absorbance spectra of "Ninghuafu" and "Donghu" samples without any preprocessing. No obvious spectral differences could be observed between the two spectra, and the spectra are highly overlapped except in the regions 940-1025 nm, 1125-1360 nm, 1670-1840 nm, and 2000-2100 nm, where "Ninghuafu" is a little higher than "Donghu", while at 1390-1435 nm as well as 2380-2480 nm, "Donghu" is a little higher than "Ninghuafu". All spectra have intense absorption bands at 1450 nm, related to the first O-H overtone, and at 1900-1950 nm, related to the combination of stretch and deformation of the O-H group in water. The small absorption band at 1690 nm might be related to the -CH3 stretch first overtone or C-H groups in aromatic compounds, at 1782 nm related to the C-H stretch first overtone, at 2266 nm likely with C-H combination bands of methanol, and at 2302 nm with combination band of C-H stretch and deformation of C-H from the -CH2 group [START_REF] Yu | Discrimination between Chinese rice wines of different geographical origins by NIRS and AAS[END_REF][START_REF] Cozzolino | Feasibility study on the use of visible and near-infrared spectroscopy together with chemometrics to discriminate between commercial white wines of different varietal origins[END_REF][START_REF] Yan | Basic and Application of Near-Infrared Spectral Analysis[END_REF][START_REF] Lu | Modern Near Infrared Spectroscopy Analytical Technology[END_REF][START_REF] Dambergs | Rapid analysis of methanol in grape-derived distillation products using near-infrared transmission spectroscopy[END_REF]. Figure 2 shows the averaged second-derivative spectra of mature vinegar from the two geographical origins, from which we can observe that the averaged second-derivative spectrum of "Ninghuafu" is higher than that of "Donghu" except for an overlap at 1900-1950 nm due to intense absorbance of water. There are obvious differences between the three spectra that are highly overlapped in figure 1, which demonstrate that features of the two spectra are enhanced after preprocessing with second derivatives.

Principal component

All the spectra of 36 "Ninghuafu" samples, 36 "Donghu" samples were preprocessed with multiplicative signal correction (MSC) and second derivatives to reduce spectral variations and enhance features of the spectra before PCA. The first ten PCs account for 99.475% of the variation in the spectra, and are used to make differentiation clearer. From figure 3, we can see that the samples were divided into two groups. The separation of "Ninghuafu" from the other was clear, but with some overlapped samples. However, the classification trends are obvious in this figure, and samples of the three brands can be separated. Figure 3 shows classification of "Ninghuafu" and "Donghu" using discriminant analysis method with obvious boundaries that include an overlapping of the two brands. It can be observed in this figure that the "Ninghuafu" and "Donghu" samples had both positive and negative scores.

From the plot, we can see that "Ninghuafu" and "Donghu" were separated perfectly with 5 overlapped samples. This result suggested that discrimination among the two brands was possible, but as the result was not ideal and PCA only represented a potential ability to discriminate the samples, discriminant analysis was employed in order to obtain a better classification.

Discriminant analysis

As a rule, the DA procedure supplies two kinds of functions: discriminant functions, which provide the separation among groups, typically in one-or two-dimensional space, and classification functions, which are used to assign a new case into one of the known groups (Nicolas et al., 2000; Yu and Ren, 1999). In this work, the former function was used to discriminate the two brands. Sixty samples (30 "Ninhuafu" and 30 "Donghu") were selected as the calibration set, and the 12 samples (6 "Ninhuafu" and 6 "Donghu") were used as the validation set. And the details were listed in Table 1. Plots of the Mahalanobis distances of all samples to "Ninhuafu" and to "Donghu" are shown in figure 4. It can be seen in these plots that the two brands were divided into two clusters. On the whole, a good result was obtained, with a classification accuracy of 95.9%. The result indicates that DA was effective for the discrimination of the two brands and produced a good classification model.

Conclusions

The result obtained in this work suggests that NIR spectroscopy, together with chemometrics methods such as DA classification models based on PCA, is a powerful technique to classify objects, and it was effective for discriminating mature vinegar samples from different different geographical origins ("Ninghuafu" from Taiyuan, and "Donghu" from Qingxu). In this case, the best classification, with accuracy up to 96.3%, was attained by DA based on sample spectra pretreated with second derivatives. This technique can acquire enough information to offer an ideal result, avoiding time-consuming, costly, and laborious chemical and sensory analysis. This result demonstrates that NIR spectroscopy has a good potential for use as an alternative technique for developing an accurate, rapid detector to discriminate mature vinegar from geographical origins. It is unknown which components account for the variation and act in the discrimination, so further study of the chemical constituents in mature vinegar should follow. At the same time, the work reported here is a feasibility study and requires further development with considerably more samples of other brands and varieties.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Average NIR spectra of 36 "Ninghuafu" and 36 "Donghu" samples without any preprocessing

Fig. 3 .

 3 Fig. 3. Two-dimensional score plot for "Ninghuafu" samples and "Donghu" sample.

Fig. 4 .

 4 Fig. 4. Classification of "Ninghuafu" and "Donghu" using discriminant analysis method

Table 1 The result of discriminant analysis method

 1

	Group	Age (Number of samples)	Number of samples Calibration set Validation set	Correctly Classified
		1-year (10)
	Donghu	3-year (10)	30	6	96.3%
		5-year (10)
		1-year (10)
	Ninghuafu	3-year (10)	30	6	95.6%
		5-year (10)

Acknowledgement

The authors gratefully acknowledge the financial support provided by Shanxi Youth Science and Technology Research Fund (No. 2009021019-3).

The paper supported by by Shanxi Youth Science and Technology Research Fund (No. 2009021019-3).


Oded Cohen

Raphael Linker
email: linkerr@tx.technion.ac.il

Amos Naor
email: amosnaor@research.haifa

Estimation of the Number of Apples in Color Images Recorded in Orchards

Keywords: artificial vision, imagine processing, yield estimation

This work presents an algorithm for estimating the number of apples on trees using images acquired with a standard color CCD camera. The proposed system is capable of correctly identifying and localizing more than 85% of the apples in the images. To achieve this high detection rate, color and texture analyses are combined together with shape analysis. In the first step, pixels with a high probability of belonging to an "apple object" are detected according to their color and texture. In the second step, "seed areas" consisting of connected sets of pixels with a high probability of belonging to an apple object are detected. Each seed area is then extended to cover the entire visible area of the apple to which it belongs. Finally, each blob is segmented into simple components that can either be combined into circles or are discarded, so that each of the resulting circles corresponds to an apple

Introduction

Early estimation of the future yield has always been a major challenge in agriculture. In orchards, such predictions are of interest to growers, packing and storage houses, and compensation funds. For growers, flowering intensity and yield forecast would be valuable at several stages of the season to optimize the intensity of chemical and hand thinning. Today, in the absence of accurate flowering and yield forecast, growers tend to perform insufficient thinning, preferring to be on the safe side rather than causing irreversibly low yield. This has significant repercussions later in the season when labor-consuming hand thinning is required to compensate for the insufficient chemical thinning. In addition, hand thinning is performed only after the process of natural thinning has ended and the fruits have reached a size where experienced growers can roughly estimated yield by visual inspection. This delay has two serious horticulturist consequences: [START_REF] Annamalai | Citrus yield mapping system using machine vision[END_REF] it results in smaller fruits at harvest, which, in addition to having a lower market value, require more time for harvesting, and (2) high crop load before the late thinning reduces flower bud initiation, which lowers the potential bloom in the subsequent season. For packing and storing houses, accurate yield forecast is needed toward the middle of the season. Such information would enable proper planning and adequate allocation of storage spaces. Finally, the compensation funds need a yield estimate early in the season when unexpected fruit drop or hail damage may occur.

Currently, yield forecasts are based on manual counting of the number of apples on selected trees. This method is extremely time-consuming and the small number of trees that can be inspected is insufficient due to the high variability of the yield that exists in apple orchards. More accurate forecast of flowering intensity and yield will increase the confidence of growers to perform early chemical thinning, which will minimize the needs for hand thinning, will increase the fruits' size and will increase potential flowering in subsequent seasons.

The present study focused on the development of a color-imaging system for automatic estimation of the number of apples present on a tree. The more general task of fruits localization on trees using artificial vision has been investigated in numerous studies, and most of these are summarized in the excellent review paper of Jimenez et al. [START_REF] Jimenez | A survey of computer vision methods for locating fruit on trees[END_REF]. These studies used either standard color cameras as in the present study (e.g. [START_REF] Annamalai | Citrus yield mapping system using machine vision[END_REF], [START_REF] Bulanon | Optimal thresholding for the automatic recognition of apple fruits[END_REF], [START_REF] Pla | Colour segmentation based on a light reflection model to locate citrus fruits for robotic harvesting[END_REF]), multispectral or hypersectral imaging (e.g. [START_REF] Kim | Apple Yield Mapping Using a Multispectral Imaging Sensor[END_REF], [START_REF] Safren | Detection of green apples in hyperspectral images of appletree foliage using machine vision[END_REF]), or thermal imaging (e.g. [START_REF] Stajnko | Estimation of number and diameter of apple fruits in an orchard during the growing season by thermal imaging[END_REF]). As noted in some of these studies, the task is especially challenging under natural illumination, and under such conditions shading and high contrast transitions are a main problem. For instance, sunlit fruits are ten times brighter than shadowed leaves, and sunlit leaves are four times brighter than shadowed fruits [START_REF] Kim | Apple Yield Mapping Using a Multispectral Imaging Sensor[END_REF]. When the fruits of interest are green as in the present study, the task is made more complex by the low color contrast between the background and the fruits.

Problem Description and Challenges

In ideal conditions, apples would appear as non-touching circular objects with smooth surface and distinctive color such as shown in Figure 1. In such a situation, finding and counting the apples within the image would be quite simple and would basically consist of finding blobs (clusters) of pixels within some predefined color range. However, such ideal situations are rarely found in images recorded in natural outdoor conditions. In such images, most apples are partially hidden by other apples and leaves, and their color is influenced by the environment, the light conditions and the photography parameters.

Fig. 1. Apples in ideal image

Light conditions have a major influence on detection feasibility. Figure 2 shows the two major effects of direct sunlight: (1) in some areas the CCD sensor is saturated and therefore any color information is lost, and (2) shading makes it impossible to recognize some of the apples as smooth surfaces. Having identified these negative factors, we tried to determine photography parameters that would eliminate or significantly reduce such problems. The best results were achieved by acquiring the images under diffusive sunlight conditions (pictures taken near sunset), and by manually underexposing the images compared to the camera's automatic setting (exposure reduction 0.7 units). Lowering the amount of light that reaches the sensor in such a manner brings the apples to the middle of the dynamic range of the sensor and avoids saturation.

Figure 3 shows a typical image obtained with such optimal photography parameters. Although to the human eye such an image may seem too dark and "flat", it is much more suitable for computer analysis than a bright and shiny image. Visual inspection of the numerous images collected during this study led to identifying the following common situations and developing an algorithm that could handle them accordingly:

Apple Clusters and Leaf Occultation

Most apples are partially hidden by other apples and/or by leaves (Figure 4)

Shading

The natural light source (sun) causes shades and non-uniform illumination. In particular, as shown in Figure 5, parts of an apple may appear much darker than the rest.

Saturation

Sunlight can be reflected more intensely from some apples or from the apple surface regions that are perpendicular to the light source and this might cause local saturation. In such regions the red, green and blue components reach their maximum value (or close to it), and all color information is lost (Figure 6). Manual under-exposure of the images eliminates most, but not all, of these cases.

Color Inference

The color of an object is influenced by the light reflected from surrounding objects. In the present case, apples that are deeper within the tree will not only appear darker but will also be more similar in color to the leaves surrounding them (Figure 7). A first set of images was taken using the fully automatic mode of operation of the camera (Fujifilm FinePix S8000fd). Eight of these images were selected as "calibration images", which were used to develop the algorithm, and another nine were selected as "validation images". About 70% of the apples present in the calibration images were marked manually to provide "apple" calibration pixels. A second set of images was taken in the following season with a different camera (Olympus C740UZ) and new photographic parameters in an attempt to overcome some of the problems identified when analyzing the first set of images:

1. The camera shutter was set manually to 0.7 units lower than the automatic camera setting in order to darken all the objects and bring the apple pixels close to the middle of the dynamic range of the sensor. 2. The pictures were taken close to sunset, when there was nearly no direct sunlight and lighting was diffusive.

Nine of these images were used to calibrate the detection algorithm, which was then tested on another eight images.

Algorithm Description

The algorithm includes three main steps which are described below. The algorithm was developed in the ImagingChef environment (a full description of this imaging development environment can be found at: www.odedcohen.com), which supported the entire development process, from reference object marking to algorithm development and result visualization

Step 1 -Apple pixels detection by color and texture

In the first step, pixels with high likelihood of belonging to an "apple object" are detected according to their color and texture. To achieve this, a K-nearest-neighbors (KNN) classifier was built using pixels belonging to "apple" and "non-apple" objects that were manually marked in the images belonging to the calibration set. Objects overlapping the apples, such as leaves, branches and the petal area were also marked in order to exclude them from the study process. Figure 8 shows some marked objects on part of a typical image. Although this is not visible in the image, the properties defined for each object include (in addition to the object type) the relative depth of the object, so that the analysis takes into account partial occultations. Each pixel in the calibration set was characterized by its Red-Green-Blue (RGB) components and its smoothness, which was estimated as the variance of the RGB components on a small area (6-by-6 pixels).

An example of the classification result is shown in Figure 9, in which the lighter pixels have a higher probability of belonging to an apple object. Two typical problems are visible in Figure 9: 1. Saturated (bright) areas get low probability, as their color is practically white (e.g. area indicated by yellow arrow); 2. Pixels that belong to the bottom side of leaves are incorrectly recognized as "most probably apple" since this side of the leaves is lighter and its color is similar to that of apples (e.g. area indicated by pink arrow). In order to minimize this problem, the camera should be located high enough and pointing downward so that only the top side of the leaves is visible.

Step 2 -Apple surface detection by seed area detection and growth

In this step connected sets (blobs) of apple pixels are detected and extended to cover the area of the apple to which they belong. Each "seed area" consists of a connected set of pixels that have a high probability of belonging to an apple object. Ideally, each apple should result in one seed area that should cover most, if not all, of the apple. In practice, parts of the apple surface might be considerably darker or lighter, because of the amount of light that reaches it and of the way this light is reflected. Such areas are misclassified as "non-apple", resulting in a smaller seed area. Also, some apples appear as split objects due to partial occultation by branches or leaves, which results in several seed areas for the same apple (Figure 10). After detecting the seed areas, these areas are expanded to contain similar neighboring pixels. This compensates for the misclassification of the very bright and very dark pixels described above. The expansion is performed by extending each blob to include neighboring pixels with low variance. Figure 12a shows a variance map of Figure 11a in which the darker pixels indicate a higher variance while smooth surfaces appear as lighter regions. Figure 12b shows the results of the seed area expansion: both saturated and dark areas around the seed areas are now included in the "apple" area. Transforming the blobs into apples requires shape analysis. For this, the contour of each blob is segmented into the following components (Figure 14):

• arcs,

• linear segments,

First set of images -images recorded in automatic exposure mode

The results are summarized in Table 1, which shows the number of correctly detected apples and the number of false positive detection (incorrect detection of apple where there is none). These results are compared to a number of apples visible to a human observer inspecting the images for about 20-30 seconds. The overall rate of correct detection is above 85% for both the calibration and validation images while the false positive rate is 12% and 19%, respectively. A detailed analysis of the results led to the following conclusions:

• As apples are much lighter than leaves, automatic setting of the exposure by the camera resulted in apple pixel values too close to the end of the dynamic range of the sensor, which often caused some of the pixels to be saturated.

•

The light was not diffusive enough. Even though the pictures were taken early in the morning, the light distribution was still not uniform enough to enable better detection.

Second set of images -images recorded with manually-lowered light exposure

Adjusting the photography parameters as detailed in Section 3 improved the pictures quality for subsequent analysis. There were no cases of color saturation and the apples appeared as more uniform. As a result, the overall rate of correct detection increased to 90% and the false positive rate decreased to less than 10% (Table 2).

Conclusions

Apple detection in images taken under natural daylight conditions has two main inherent difficulties: (1) the natural light is not diffusive enough and might cause shades and saturation, and (2) the apples have different shapes and colors, overlap other objects, and are rarely fully visible. Nonetheless, the present algorithm is capable of detecting correctly more than 85% of the apples present in an image obtained using the automatic settings of a standard color camera. The performance can be improved by taking care of recording the images under diffuse light conditions and manually lowering the exposure in order to avoid light saturation. In such cases close to 90% of the apples are correctly detected. The proposed algorithm should be further validated using a larger dataset and its extension to other apple varieties should be considered.

Fig. 2 .

 2 Fig. 2. Apples exposed to direct sunlight that causes saturation of the CCD and strong shadows

Fig. 3 .

 3 Fig. 3. Optimal image recorded under diffuse light (close to sunset) and with manually reduced exposure

Fig. 4 .

 4 Fig. 4. Detail of a typical image with partial occultation of some of the apples

Fig. 5 .

 5 Fig. 5. Detail of an image showing non-uniform light distribution on the apples

Fig. 6 .

 6 Fig. 6. Detail of an image with saturated regions in which all color information is lost

Fig. 7 .

 7 Fig. 7. Detail of an image showing the large color differences that exist between the apples on the outside of the tree and those deeper within the tree

Fig. 8 .

 8 Fig.8. Detail of an image with objects marked using the ImagingChef software. In addition to the apples (indicated by circles), regions that should be ignored at the calibration stage (such as apples that are barely visible or parts of apples that are occulted by leaves and branches) are marked by rectangles or ellipses.

Fig. 9 .

 9 Fig. 9. (Left) Original image and (Right) image showing the classification results, in which the lighter pixels have a higher probability to belong to an apple. The yellow and pink arrows point to typical incorrect classifications as explained in the text

Fig. 10 .

 10 Fig. 10. Detail of an image that shows that a single apple may contain more than one seed region (purple contour) as a result of partial occultation Figure 11 shows typical results of the seed areas detection procedure. It can be seen that both saturated and darker pixels of the apple surfaces are not included in the seed areas.

Fig. 11 .

 11 Fig. 11. (Left) Typical image. (Right) Detected "seed areas" (in orange)

Fig. 12 . 4 . 3 Step 3 -

 12433 Fig. 12. (Left) Variance map of Picture 11a. (Right) Seed areas obtained by expanding the areas shown in Picture 11b using the variance map shown in Picture 12a

Fig. 13 .

 13 Fig. 13. Typical seed areas

  In the next step, the arcs are grouped into circles and at the end of this stage each circle indicates one apple.

Fig. 14 .

 14 Fig. 14. Contours obtained after segmenting the contour of each seed region into arcs, linear segments and anamorphic segments

Table 1 :

 1 Results of the analysis of the images recorded using automatic light exposure

		Calibration images
		Number of	Number of	Number of
		apples	apples	false
		visible	detected	detections
	Image 1	31	26	2
	Image 2	29	25	2
	Image 3	45	34	1
	Image 4	37	33	8
	Image 5	22	17	4
	Image 6	31	26	3
	Image 7	19	19	2
	Image 8	49	44	9
	Total	263	224 (85%)	31 (12%)
			Validation images
	Image 1	44	44	6
	Image 2	26	23	4
	Image 3	40	38	12
	Image 4	19	17	9
	Image 5	22	21	5
	Image 6	60	53	4
	Image 7	47	37	10
	Image 8	38	33	0
	Image 9	44	34	11
	Total	296	256 (86%)	55 (19%)

Table 2 :

 2 Results of the analysis of the images recorded under diffuse light and after setting manually the camera shutter 0.7 units lower than prescribed by the automatic setting

		Calibration images
		Number of	Number of	Number of
		apples	apples	false
		visible	detected	detections
	Image 1	66	61	4
	Image 2	68	61	7
	Image 3	29	28	7
	Image 4	28	25	5
	Image 5	37	33	2
	Image 6	13	13	3
	Image 7	36	33	3
	Image 8	58	53	7
	Image 9	29	28	1
	Total	364	335 (92%)	39 (11%)
			Validation images
	Image 1	57	47	12
	Image 2	48	39	2
	Image 3	60	52	4
	Image 4	17	14	2
	Image 5	31	29	2
	Image 6	24	17	1
	Image 7	53	48	1
	Image 8	62	57	7
	Total	295	256 (87%)	19 (6%)


Xueliang Fu

Changyou Li

Buyu Wang

Honghui Li

Hailei Ma

Dongnan Zhu

The Monitoring System of Water Environment Based on Overlay Network Technology 1

Keywords: :WSN, Overlay network, ETL, Multidimensional data cube

At present, although methods of automatic, periodical manual data collection have been adopted in some areas, there are many problems existing in these methods. This paper proposes a new framework for online monitoring of water environment based on overlay network technology. Especially, we build the multi-level overlay network, which consists of the GPRS networks, mobile networks and internet networks. Following this, a multi-dimensional data cube for water environment is established using the ETL process with the input of complex heterogeneous data collected. Thereby, the framework of data center is established for on-line early warning of water environment and data analysis and processing. This framework has been put into use in Wuliangsuhai for the on-line real-time water monitoring tests. The results show that the overlapping network architecture is effective and on-line analysis and early warning of heterogeneous data is efficient.

Introduction

Currently, water-quality monitoring in Wuliangsuhai is mainly performed manually [START_REF] Seelig | Relations of remote sensing leaf water indices to leaf water thickness in cowpea, bean, and sugarbeet plants[END_REF]. Researchers in related institutes reside in Wuliangsuhai, and gather water-quality data manually during the ice-absent period. The main method to gather data relies on portable instruments. It is difficult to meet the needs of comprehensive monitoring, analysis and early warning in Wuliangsuhai for water-quality. Therefore, the comprehensive development of research and eco-social benefit is seriously restricted [START_REF] Ross | Instrument with integrated sensors for a rapid determination of inorganic ions[END_REF].

There are mainly four issues in the above traditional method. Firstly, researchers had to be outside for long time, thus, they cannot conduct large-scale experimental analysis using tools and computing environment. Secondly, the continuity of data collection is poor, and a sample space is very limited. In the ice-absent period, one can only collect data once or twice per day (once at noon, or twice in the morning and evening respectively). One cannot collect data at night as well as in freeze-up period. So, the number of the achieved water-quality data is no more than 200 records. Agriculture Therefore, it is impossible of widely data mining and analysis owing to the significant shortage of data. Moreover, the quality of the collected data cannot be guaranteed. Thirdly, using the manual-collecting method, the quality of data is effected strongly by the professional proficiency of the researchers and the usage of the equipment. As a result, dirty data may appear and cannot be analyzed and rectified using other timesequential data. Finally, the collection area is very limited. Due to the broad water area of the Wuliangsuhai Lake, it is impossible for the residing researchers to collect data at the same time from multiple locations far away from each other. As a result, there is no way to analyze and compare data in terms of time scale and location.

The current trends to resolve these issues is to use the framework of cloud services, with which data is online automatically assembled, real-time transmitted using overlay networks. Decision of knowledge is made in a data center.

Wireless sensor networks (WSN) are the task-oriented wireless networks which consist of a number of wireless sensor nodes [START_REF] Tatyana | A sensitive nitrate ion-selective electrode from a pencil lead[END_REF]. WSN integrates multiple area technology, including sensor technology, embedded computing technology, modern networking, wireless communication technology, distributed information processing technology and others. In WSN, various micro sensors take responsibility of on-line monitoring target; the embedded computing resources take care of processing data obtained by sensors; the related information is send to the remote user data center using wireless communication networks. This technology can be broadly applied in the military defense, industrial and agricultural controls, urban management, biological medicine, environmental monitoring, disaster relief, antiterrorism and remote control in dangerous areas. It is attractive in both academy and industry [START_REF] Wang | A clustering Algorihm Based on Latent Semantic Model[END_REF].

This paper studied the design of monitoring system for environment of water based on the overlapping networks in order to deal with the issues about Wuliangsuhai water-quality monitoring.

Main Contributions

Architecture

Our proposed system is mainly composed of overlapping network communication subsystem, the water-quality sensor acquisition subsystem, multi-dimension data analysis subsystem [START_REF] Cooley | Remote Sensing of the Coastal Zone of Tropical Lakes Using Synthetic Aperture Radar and Optical Data[END_REF] and solar power subsystem. The issues are addressed using multi-layer overlay network technologies on inter communication among networks, such as 485 network, GPRS network, PSTN network and computer networks. In this system, the task-oriented data-collecting network is constructed using water quality sensors and data acquisition equipment based on wireless sensor network technology; the multi-dimensional data analysis center are established based on OWB and ETL technologies for early warning and comprehensive analysis of instantaneous data and historical data; solar power system is employed for 24-hour uninterrupted power supply all day long and thus provide guaranteed environments for the good performance of the entire system. System architecture is shown in Figure 1.

The Monitoring System of Water Environment Based on Overlay Network Technology 3

Overlay Network Design

A networking framework which enables multi-network integration must be explored in order to realize a distributed resource management platform. This platform can provide for the application layer reliable and efficient resource searching and localization service through monitoring link status, shielding network failures and changes. Also, this platform can provide for application service such as network routing optimization through node detection of network paths and performance.

The real-time signal data on water quality parameters is achieved through the data acquisition subsystem. The collected signal data is first send to a GPRS-wireless hybrid communication system via a RS485 network, and then transfer to the data center with fixed IP address by the GPRS wireless network. In the data center, the signal data is translated, analyzed, processed and stored. The TCP multicast protocol, i.e., "center to multi-point", is adopted for reliable and transparent data transmission.

Multiple-layer hybrid overlay networks are constructed between each sensor node, which registers as an online communication node and the super node in the data center with a fixed IP. The super node manages the whole WSN. Each facility, which is eligible for the WSN, can log in, exit and abnormal exit through this subsystem. This subsystem is important for the establishment of a reliable network, which provides bidirectional data transmission channels and communication links for the transmission of signal data and control signaling.

The layered hybrid network architecture is exploited in the whole network framework. The whole network is divided into four layers: the service layer, the core layer, the access layer and the link layer.

Data Stream Processing

The Data center architecture is shown in Figure 2. The data stream processing includes two steps: receiving data stream and sending data stream, which are described as follows:

1) Receiving data stream i. The on-line real-time data on water quality parameters is send to the data center in the sampling frequency using the data acquisition subsystem. The acquisition subsystem consists of water-quality-parameter sensors, filters and A / D converters.

ii. After receiving the real-time data, the data center analyzes the data logically, eliminate the dirty data, and store the qualified data in the transient database.

iii. Water-quality analysis-type database or data warehouse, which is suitable for analysis and statistics, is built from the data in transient database using the ETL tools combined with the functionalities of operations / scheduling of enterprise database.

iv. Various analysis products needed by users can be generated using the database products for analysis (BI data warehouse reporting engine) together with the waterquality analysis-type database or data warehouse.

2) Sending data stream i. A user can log in the data center using the mobile terminal equipment (e.g., mobile phones, PDA, etc.) or the water quality data viewer on a fixed terminal.

The Monitoring System of Water Environment Based on Overlay Network Technology 5

Through the data center, the user can transmit the acquisition subsystem in order to set the sampling period and other parameters of the acquisition system.

ii. With the help of the data center for sending the control information, an eligible user can turn on or off the power supply system, or remote manage the power supply system.

Experiment setup and discussion

Firstly this paper deployed the PH-value and ORP and oxygen-content sensors in the Wuliangsuhai Lake. Then, the real-time, automatic data acquired with adjustable sampling frequency has been achieved using the proposed water-quality sensor acquisition subsystem, the overlay network transmission subsystem, and the data storage and processing subsystem. Thereby, the water-quality parameter data can be processed and analyzed in a very short time.

The software in PLC (Programmable Logic Controller) exploits the means of center-to-multi-point communication in the water-quality-parameter signal acquisition subsystem in Wuliangsuhai; the Baud Rate is 9600; the survivable mechanism is achieved through the heart-rate packet; the communications signaling employs ASCII-code signaling. The connections between the PC software and the data center are established using TCP/IP protocol. The unity of the user programming interface is approached through encapsulating the Socket package in the system level and the integrated converged communications from the data link layer, the network layer and the application layer.

The data center is built based on Oracle database, where the transaction concurrency mechanism and trigger mechanism are used, and the data OLTP is achieved through the job scheduling. Oracle DWB is used to build the data warehouse which provides analysis environments. All analysis products and user UI are implemented using the B/S framework; the system is built using the SSH enterprise multi-layer computing framework.

After the deployment in practice and the simulations, it is found that the proposed overlay network framework is sensitive very much to the strength of the data signals.

The instantaneous packet loss may suffer when the signals become weak, which results in dirty data. However, due to the relatively high sampling frequency, this can be ignored in the context of the research and the application on water environment monitoring only if the data is accurate in the time granularity of an hour.

Conclusion

This paper proposed and implemented the Monitoring system of water environment based on overlay network technology in order to deal with the issues on the current water environment monitoring in the Wuliangsuhai Lake. Through the field deployment of our proposed system, simulation results verify that the effectiveness of our design. It solves a series of issues existing in the current monitoring method mentioned before. The future work is to investigate the issues

The research is supported by Chinese Natural Science Foundations (50969005,40901262) and by Specialized Research fund of High Education for Inner Mongolia (Njzy08046)

The Fourth International Conference on Computer and Computing Technologies in Agriculture

Agriculture

related to data processing and data reliable transmission both in academy and in practice.


Gailian Li

, Tingxiang

Chunxia Hu

Analysis of the heat transfer performance of vaporcondenser during vacuum cooling

Keywords: Heat transfer performance, Vapor-condenser, Cold trap, Vacuum cooling

The heat transfer performance of vapor-condenser has been studied under different temperature of cold trap and different thickness of the frost layer in this paper. The relationship between dimensionless number and is obtained. The results show that the capturing efficiency of cold trap increases with the decrease of the surface temperature of vapor-condenser. The frost accumulated on the surface of vapor-condenser can cause the overall heat transfer coefficient decrease, which has a negative effect on heat transfer of vaporcondenser.

has an effect on the heat transfer of vapor-condenser during vacuum cooling.

Introduction

Vacuum cooling is a rapid evaporative cooling method. Vacuum cooling has been successfully used to cool vegetables and flowers since the 1950s [1]. In the recent years, for the safety of foods, a rapid cooling treatment after cooking process should be used to minimize the growth of surviving organisms. Compared with the conventional cooling methods, vacuum cooling has many advantages. Therefore, many researches have highlighted the applications of vacuum cooling for the cooked meats [2][3][START_REF] Desmond | Effect of rapid and conventional cooling methods on the quality of cooked ham joints[END_REF]. In addition, heat and mass transfer characteristics during vacuum cooling have been investigated. Predictive models can provide much valuable information for the cooling process of large cooked meat joints under broad experimental conditions within a short time. Wang and Sun have developed a mathematical model for describing the vacuum cooling process of the large cooked meat joints [5][START_REF] Wang | Modelling vacuum cooling process of cooked meat-part 2: mass and heat transfer of cooked meat under vacuum pressure[END_REF][7].

A vacuum cooler is a machine to maintain the defined vacuum pressure in a sealed chamber, where the boiling of the water in the cooked meats occurs to produce the cooling effect. Theoretically, only the speed of vacuum pump is high enough to produce the defined vacuum pressure in the vacuum chamber. However, at a low

Nomenclature

D -Gas constant, ; 1 1 - -⋅ ⋅ K mol J - λ Thermal conductivity, ; 1 1 1 - - - ⋅ ⋅ ⋅ s K m J - t Time, ; s - σ Stefan's constant, ; 4 2 8 10 7 . 5 - - - ⋅ ⋅ × K m W - A Area, ; 2 m - ε Emissivity; - d Diameter, ;; m λ -Mean free path of free molecular, ; m - α Heat transfer coefficient ; Subscripts 1 2 - -⋅ ⋅ K m W - i o h h , Enthalpy, ; -Frost layer; 1 - ⋅ kg J fr - p C Specific heat at constant pressure, ; v -Vapor; 1 1 - -⋅ ⋅ K kg J - t Time, ; s - i Inlet; - t F Thermal accommodation factor; - o Outlet; - F Tangential momentum accommodation factor; - c Coolant;
pressure, the volume ratio of steam and water is very large. For example, when the pressure is 1073 Pa, the corresponding saturation temperature is 8℃, the specific volume is 120.851 kg m 3 . If the entire vapor is evacuated only through the vacuum pump, the speed of vacuum pump should be very large, many vacuum pumps are required in the vacuum cooler, which is obviously unsuitable. In order to remove the large amount of water vapor and keep the cooling cycle within a reasonable length of time, the vapor-condenser is used to economically and practically handle the large volume of water vapor by condensing the vapor back to water and then draining it through the drain valve. The vacuum pump and the vapor-condenser in the vacuum cooling system are used to remove the water vapor evaporated from the cooked meats. Generally, the temperature of the vapor-condenser is about -30℃~-50℃ during vacuum cooling. The large temperature difference exits between the surface of vaporcondenser and the water vapor in the vacuum chamber. Consequently, the water vapor will become the frost at the surface of vapor-condenser. The frost formation on the cold surface below 0℃ acts not only as a thermal insulator between the surface and the water vapor, but also significantly reduces the heat transfer performance of vaporcondenser, which can result in the decrease of the capability of capturing water vapor in the vapor-condenser. In order to improve the capturing efficiency of the vaporcondenser, the evaporation temperature of refrigerant in the vacuum cooler must be lowered. However, the lower evaporation temperature will add the cost of vacuum cooler and energy consumption. Therefore, it is very important to investigate the heat transfer performance of vapor-condenser for designing and optimizing the vacuum cooling system. The phase change heat transfer theory under vacuum pressure is hardly studied. Hong and Leena [START_REF] Chen | Modeling frost characteristics on heat exchanger fins: Parts I, Numerical Model[END_REF] have modeled the frost characteristics under atmosphere pressure. In the current study, the heat transfer performance of vaporcondenser in the vacuum cooler is investigated. Moreover, the factors of affecting the heat transfer performance of vapor-condenser are also analyzed.

Experimental apparatus

The laboratory-scale vacuum cooler as shown in Fig. 1 was built by Shanghai Pudong Freezing Dryer Instruments Co. Ltd. (Shanghai, China). Vacuum cooler has four basic components: a vacuum chamber, a vacuum pump, a vapor-condenser and a refrigeration system. The vapor-condenser is an evaporator in the refrigeration system and a condenser capturing water vapor evaporated from the cooked meats during vacuum cooling. The cooling coil of vapor-condenser is set up in a stainless cylindrical steel, which is enclosed with 30 mm thickness polyurethane foam to prevent heat transfer. The stainless cylindrical steel with vapor-condenser is defined as cold trap. The water vapor is evaporated from the cooked meats during vacuum cooling. The vapor-condenser and vacuum pump removes the water vapor and air to reduce the pressure in the vacuum chamber. Because of the large temperature difference between the cold trap and the water vapor, the large amount of water vapor can enter into the cold trap. One part of water vapor is condensed into water by liquefaction, and the other part of water vapor become frost on the surface of vaporcondenser by solidification. A set of T-type copper-constantan thermocouples with an accuracy of 0.1 ± ℃ are used to record the temperature of the cold trap. The mass flux of coolant is measured through supersonic flowmeter (UFLO2000P, USA). The vacuum pressure is measured through the pressure transducer (CPCA-130Z) with an accuracy of Pa, the range of pressure transducer is 10 Pa~10 Kpa. The mass of the sample and water are measured through the electric balance (JA12002, made in China). Thermal conductivity of frost layer is measured by an unsteady state method using a line heat thermal conductivity probe, based on the design of Sweat as described by Scully [13,14]. The thickness of frost layer is directly determined by a micrometer having a 0.1mm resolution. 0.5 ± 1-bleeding valve; 2-weight sensor; 3-sample; 4-thermal couple; 5-pressure sensor; 6-vacuum chamber; 7-electronic balance; 8-compute; 9-temperature controller; 10-coolant outlet; 11coolant inlet; 12-cold trap; 13-vacuum pump; 14-pressure controller; 15-I-7018P module The formation mechanism of the frost on the surface of vapor-condenser Fig. 2 shows the formation process of frost. The sensible heat is transferred from the water vapor in the cold trap to the frost surface by the temperature difference driving force between the water vapor and the frost surface. Some of the transferred moisture deposits on the frost layer, causing the frost layer to grow. The remainder diffuses into the frost layer. The heat of sublimation caused by the phase change of the added frost layer is transferred through the frost layer. The latent heat and sensible heat transferred from the water vapor are then transferred through the frost layer by conduction. The water vapor diffusing into the frost layer changes phase within the frost layer. The frost density increases as a result of this process.

The frost layer is a porous medium composed of ice crystal and air. The ice crystal has different shapes during the formation of the frost layer. Ice crystal shapes are classified into main forms: plate-like forms and column-like forms. The microscopic structure of ice crystal is shown as in Fig. 3 [START_REF] Webb | 13. ds.), Engineering 14. ges in beef burgers and mashed[END_REF].

[ ] ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - + = ice r f ice v v vs v ice r f ice v fr fr vs fr v v P h D v v RT h Q m ρ ρ ρ ρ λ 1 ) ( 1 ) ( 5 . 0 2 0 & (1)
Where is the refrigeration load of vapor-condenser; 0 Q vs h is the sublimation heat of ice; R is the gas constant; fr T is the surface temperature of the frost layer;

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = dr d dt d r r m dt d r f r f r f ρ ρ 4 2 & (3)

Heat transfer in cold trap

The temperature of cold trap is about -30℃~ -60℃. The water vapor evaporated from the cooked meats will become the frost at the surface of vapor-condenser in the cold trap. The frost layer at the surface of vapor-condenser gets thicker and thicker with the increment of time. The frost layer is a porous structure composed of ice crystal and air pores. Moreover, the porous structure contains a low thermal conductivity, which reduces the thermal conductivity of the frost layer. Finally, the frost layer results in a significant heat transfer resistance from the water vapor to the surface of vapor-condenser in the cold trap. The relationship between the heat flux and the thickness of the frost layer can be expressed:

δ λ T q r f Δ - = (4)
Where fr λ is the thermal conductivity of the frost layer; δ is the thickness of the frost layer.

Heat transfer coefficient, is an important index to evaluate the heat transfer performance. Heat transfer coefficient of cold trap is measured by the qusai-stable method in this experiment. The total heat transfer coefficient is expressed as follows:

k m T A Q k Δ ⋅ = 0 (5)
Where A is the total area of heat transfer; the logarithmic temperature difference, , can be expressed as:

m T Δ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - = Δ e o e i o i m T T T T T T T ln (6)
Where is the evaporation temperature.

e T

On the other hand, the total heat transfer coefficient is also theoretically determined by:

2 1 1 1 1 d d k c fr v ⋅ + + = α λ δ α (7)
Where v α is the heat transfer coefficient of vapor in cold trap; c α the heat transfer coefficient of coolant in coil of vapor-condenser; and are respectively inner and outer diameters. The refrigeration load of vapor-condenser, , in Eq. ( 1) and Eq. ( 5) can be calculated by the enthalpy difference of refrigerant between inlet and outlet.

0 Q ( ) i o c h h m Q - = & 0 (8)

Where

is the mass flux of refrigerant; , is respectively the enthalpy of refrigerant in outlet and inlet.

c m & o h i h 0
Q is also calculated through heat transfer of gas in cold trap. The expression can be given by:

( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + - + - = ∫ ∫ ∞ ∞ ∞ ml fr ml T T T T vs v p r pf r f fr fr v h dT C dT C m T T A T T A Q & 4 4 0 σ ε α (9)
Where , is the specific heat of the frost and water vapor; is the gas temperature in cold trap; is the sublimation temperature.

r pf C v p C ∞ T ml T
The vacuum pump and vapor-condenser can cause the reduction of pressure in cold trap. When the gases in cold trap are at low pressure, the slip flow occurs. The relative importance of effects due to the rarefaction of a gas in cold trap can be indicated by Knudsen number ( ), a ration of the magnitude of the mean free molecular path ( Kn λ ) in the gas to the characteristic dimension ( ) in the flow field. When the slip flow occurs in the cold trap, the gas adjacent to the surface no longer reaches the velocity or temperature of the surface. The gas at the surface of the frost has a tangential velocity and it slips along the surface. The temperature of the gas at the surface of the frost is finitely different from the surface temperature of the frost layer, and there is a jump in temperature between the surface of the frost layer and the adjacent gas. The energy and momentum equations in cylindrical coordinates can be written as [12]:

L ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ = ∂ ∂ r T r r r C x T u p 1 ρ λ (10) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ = ∂ ∂ r u r r r x P 1 μ (11)
The slip velocity as a function of the velocity gradient near the wall of cold trap can be expressed as:

0 0 2 r r r r r u F F u = = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ - ⋅ = λ (12)
The temperature jump in slip flow at the wall of cold trap can be written as:

0 Pr 1 2 2 r r t t w v r T F F T T = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ⋅ + ⋅ - = - λ γ γ (13)
Where λ is mean free path of water vapor; F is the tangential momentum accommodation factor;

is the thermal accommodation factor; t F γ is the ratio of specific heat; Pr is the Prandtl number;

is the radius of cold trap.

r

The Nusselt number can be given from Eq. ( 10) and ( 11):

0 2 1 Pr 1 2 2 24 6 1 6 6 1 36 11 48 r F F Kn Kn Kn Kn Nu t t ⋅ ⋅ + ⋅ - ⋅ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + + - = λ γ γ (14)

4

Result and discussion

The influence of the frost layer on heat transfer in cold trap

Fig. 4 shows the experimental data for the thermal conductivity of the frost. The time ranges, for which the data were taken, are given on Fig. 4. When the experimental temperature of cold trap is -45℃, the average thermal conductivity of the frost layer during vacuum cooling is 0.1072 . It can be also found from Fig. 4 that the thermal conductivity of the frost layer increases with the increment of time. This is because the density and the thickness of the frost layer depend on the temperature of the frost and the time. The lower temperature, the more water vapor in cold trap becomes the frost on the surface of vapor-condenser and diffuses into the frost layer, which can increase the density and the thickness of the frost layer. The relationship between the thermal conductivity of the frost layer and the density of the frost layer has been shown in Eq. ( 2). Fig. 6 shows the capturing efficiency of cold tw p. It can be found that the lower temperature of cold trap is, the higher the capturing efficiency of cold trap is. When the temperature of cold trap is about -55℃, the capturing efficiency of cold trap is above 90%. However, if the temperature of cold trap decreases to -40℃, the capturing efficiency of cold trap is about 20%. Therefore, the temperature of cold trap should be very low so that the cold trap can capture more water vapor. On the other hand, if the temperature of cold trap is too low, some water vapors become the frost on the surface of vapor-condenser and on the wall of cold trap, the thick frost layer has the low thermal conductivity. Which increases the heat transfer resistance between the surface of vapor-condenser and the water vapor. With the increment of thickness and density of the frost layer, the heat transfer resistance between the surface of vapor-condenser and the water vapor increases, the temperature of cold trap become high so that the capturing efficiency of cold trap gets more and more low. g. 7 gives the variation between Nu and Kn in the different thermal accommodation factor. Kn is correlated with t acuu pressure. It can be found that Nu decreases when the vacuum pressure in ld trap decreases. This is because convection hea n reduce at low vacuum pressure. Nu at the high thermal accommodation factor (1.0) is higher than that at the low thermal accommodation factor (0.8). The variation of Nu is opposite to the variation of Kn . he v m co t ca

The heat transfer pe temperature of cold trap e temperature of cold trap is, the more water vapor is captured. At the same time, because of the low temperature in cold trap, water vapor become frost at the surface of vapor-condenser and on the wall of cold trap. The frost layer is a porous medium composed of ice crystal and air. The low thermal conductivity of the frost layer has a negative effect on heat transfer of vapor-condenser. When the accumulated frost at the surface of vapor-condenser becomes thick, the capturing efficiency of cold trap will decrease. In addition, the relationship between dimensionless number Nu and Kn is obtained, the variation of Nu is opposite to the variation of Kn .

KNO LEDGEMENTS

Education (P. R. China).

1. Briley, GC. 53 (2004).

a rapid vacuum cooled cooked beef product. Journal of Food Engineering . 47, 139~147 (2001).

3.

Burfoot, D., Self, K. P., Hudson, W. R., Wilkins, T. J., & James, S. J. Effect of cooking and cooling method on the processing times, mass losses and bacterial condition of large meat join

γ

  The ratio of specific heat;

Fig. 1 .

 1 Fig. 1. Schematic diagram of the vacuum cooler system

Fig. 2 .Fig. 3 .

 23 Fig. 2. The formation process of the frost

  ρ are respectively density of frost layer and ice; v P is the partial pressure of water vapor; v D is the diffusivity of water vapor; fr λ is the thermal conductivity of the frost layer, the expression is as follows [frost can change during the formation of frost.

Fig. 4 .

 4 Fig. 4. Thermal conductivity of the frost layer

Fig. 5 .

 5 Fig. 5. The influence of thickness of the frost layer on the heat transfer co

Fig. 6 .

 6 Fig. 6. The influence of temperature of cold trap on the capture efficiency

Fig. 7 .

 7 Fig. 7. The relationship between Nu and Kn