Lorenz Brenner

Frank Tillenkamp
email: frank.tillenkamp@zhaw.ch

Christian Ghiaus

Refrigeration machine modeling for exergy-based performance and optimization potential evaluation of chillers in real field plants Modélisation de machines frigorifiques pour l'évaluation exergétique de performance et du potentiel d'optimisation des refroidisseurs dans des installations réelles

Keywords: optimization potential assessment method, modeling, refrigeration plant, refrigeration machine

To investigate and optimize a refrigeration system, the behavior at various operating conditions must be known or determined. The performance and improvement possibilities may then be inferred from measurement data and compared with corresponding performance key figures. These values are typically referred to normal conditions and it is usually unknown which ones represent an adequate operation. However, it is relevant for refrigeration plant operators to have reference values for a large range of operation conditions as a baseline for determining the obtainable improvements. The present work proposes the application of refrigeration machine models for increasing the range of applicability of the exergy-based optimization potential index method. Four different modeling approaches are evaluated and discussed: equation-fit, physical lumped parameter, refrigeration cycle and artificial neural network based models. By using experimental data of a real field plant, the artificial neural network approach revealed the best performance for the present application. The practical usage of the improved evaluation method is shown for the subsystem refrigeration machine on a real field installation as a case study. With the introduced additional limits, the interpretability of the results is increased.

The distinction between adequate (technical requirements exceeded), acceptable (technical requirements fulfilled) and inadequate (potential for improvement) operation according to the state of the art in technology is straightforward, which is important in practice.

The model implementation in the optimization potential index (OPI) is shown.

Additional technical limiting values for the OPI are proposed.

The performance and optimization potential of refrigeration machines is revealed.

The application of the method is shown on a real field plant as a case study.

Nomenclature

Introduction

A significant amount of the energy consumption in buildings is due to the heating, ventilation and airconditioning (HVAC) systems. Depending on the climatic region, about 30 -70% of the total energy consumption in air-conditioned buildings is due to systems for thermal comfort applications [START_REF] Vakiloroaya | A review of different strategies for HVAC energy saving[END_REF]. Thereby, more than 50% of this energy is consumed by refrigeration machines and plants [START_REF] Wang | A practical approach to chiller plants' optimisation[END_REF]. Consequently, the optimization of such systems may significantly reduce the power consumption in air-conditioned buildings.

In order to investigate and optimize a refrigeration system, the behavior at various operating conditions must be known or determined. Exergy analysis is seen as an appropriate approach, as the electrical energy consumption of the refrigeration machines and auxiliary devices is the main expense in refrigeration plant operation. This represents high quality energy which is degraded along the processes in the whole system.

In refrigeration plants the interest extends from the refrigeration machine itself to the neighboring hydraulic circuits, since a well-suited hydraulic integration of the chiller and an adequate operation of auxiliary devices is crucial in order to achieve a high efficiency of the plant. Furthermore, it is relevant for refrigeration plant operators to have reference values for determining the obtainable improvements. Interpretation of the results of the procedure should not require specialists to evaluate if there is a need to take action or not. For this reason, a practice-oriented evaluation method for refrigeration plants in air-conditioning applications based on exergy analysis and technical standards as baseline was introduced [START_REF] Brenner | Optimization potential index (OPI): An evaluation method for performance assessment and optimization potential of chillers in HVAC plants[END_REF][START_REF] Brenner | Exergy performance and optimization potential of refrigeration plants in free cooling operation[END_REF]. The method allows assessing the performance and the optimization potential with respect to the state of the art in technology and can be widely applicable in practice with most common measuring equipment of real field plants. With regard to this, the refrigeration plant is split into reasonable subsystems. This allows an individual assessment and leads to a reduction of the required measurement variables but still ensuring a sufficient level of detail. The optimization potential index (OPI) [START_REF] Brenner | Optimization potential index (OPI): An evaluation method for performance assessment and optimization potential of chillers in HVAC plants[END_REF] reveals the potential for improvement with respect to the technological baseline at a glance regardless the complexity of the system (see section 2 for details).

While for most subsystems enough data is available to derive the technological baseline from standards, it is difficult for the situation of refrigeration machines in real field plants. No specific baseline values for the electrical power consumption of the compressor are available in technical standards because it depends on a variety of variables. Consequently, further research is required to specify additional reference values for the baseline in order to apply the OPI method in a wider range of operation conditions. Ideally, these values are determined with comprehensive and representative measurements. However, installing measurement equipment is always bound to investment costs, and therefore, most of the refrigeration machines are instrumented with the indispensable sensors necessary for the operation. These variables are rarely logged for monitoring purposes, as the secondary side temperatures and cooling load is of main interest to the plant operators, e.g.

for cost allocation of the connected consumers. The latter is one reason why much attention was paid to refrigeration machine models, as it is a suitable method to determine missing quantities or to predict the system behavior with obtainable data.

s u b m i s s i o n p r e p r i n t

Consequently, various studies handle the topic of refrigeration system analysis and optimization with the aid of numerical models of refrigeration machines or other field plant devices. Among others, Shan et al.

proposed an improved chiller sequence control strategy for refrigeration plants with centrifugal chillers applying a multi-linear regression model [START_REF] Shan | Development and validation of an effective and robust chiller sequence control strategy using data-driven models[END_REF]. Typically, an optimal sequencing of the refrigeration machine operation can improve the plant efficiency and, therefore, reduce the electrical energy consumption. The authors showed an energy saving potential of 3% in comparison to the original control strategy. Wei et al. investigated a chiller plant with four refrigeration machines, four cooling towers and two cold water storages [START_REF] Wei | Modeling and optimization of a chiller plant[END_REF]. A datadriven approach was chosen to model the plant and subsequently, to ameliorate the operating conditions. The model was applied with measurement data of two days, where the authors demonstrate an energy consumption reduction of approximately 14%. Wang et al. handled the topic of chiller plant optimization to reduce the power consumption of the plant, while applying artificial neural network models to simulate the refrigeration machines and the cooling towers [START_REF] Wang | A practical approach to chiller plants' optimisation[END_REF]. Another study focused on the energy optimization of a multi-chiller plant with cooling towers in a multistory office building [START_REF] Thangavelu | Energy optimization methodology of multi-chiller plant in commercial buildings[END_REF]. Energy models for the chillers, cooling towers and auxiliary devices were developed, in order to improve the energy utilization. For the optimization routine, the building load and ambient air conditions were used as inputs. The optimization delivered ideal on / off strategies of the equipment as well as chilled and cooling water conditions. By carrying out three case studies, the authors identified an average energy saving of 20% for small chiller plants and up to 40% for moderate sized refrigeration systems.

With the stated need for research, the present work proposes the application of refrigeration machine models to introduce additional technical baseline values for the optimization potential index in the subsystem refrigeration machine. Four different modeling approaches are investigated and discussed: equation-fit, physical lumped parameter, refrigeration cycle and artificial neural network based models. All models are generally applicable with commonly measured quantities in field plants as input variables in order to ensure the practical application of the evaluation method. Their usage and performance is evaluated with measurements from a real field installation, where the best suiting model is then applied to calculate the optimization potential index. As a cases study, the application of the assessment method for the subsystem refrigeration machine is exemplified with experimental data gained from a real field plant.

Assessment method

General optimization potential index definition

The optimization potential index is given by [START_REF] Brenner | Optimization potential index (OPI): An evaluation method for performance assessment and optimization potential of chillers in HVAC plants[END_REF]:

OP I = 1 - 24 h t=0 h B * in 24 h t=0 h B in (1)
with B in and B * in the actual and reference exergy input. The optimization potential index relates the actual input of each subsystem with a reference one according to the state of the art in technology, while the same s u b m i s s i o n p r e p r i n t output is achieved. Therefore, the key figure indicates how the real system would behave in comparison to a reference system in exactly the same situation. As an approach, these reference values are derived from technical standards, as they are usually specified in tenders or contracts and should be fulfilled at the stage of commissioning. Also, these values represent an achievable technological baseline, which can depend on the requirements in different countries or regions. As the cooling load typically follows a daily rhythm, the key figure is evaluated on a daily basis, e.g. for a daily check.

The interpretation of the results is straight-forward also for non-specialists, which is important in practice:

If the actual effort is larger than the reference, an optimization potential is present which is indicated with an OPI greater than zero. Conversely, an OPI lower than zero indicates an adequate operation of the system where the technical requirements are exceeded. However, it is unknown which values close to OPI = 0 (technical requirements fulfilled) still represent a permissible operation, since the boundary between the two operating states is sharp (see Fig. 1a). By introducing an additional acceptable limit OP I acc , it can then be distinguished between adequate, acceptable and inadequate operation (see Fig. 1b), which yields a better interpretability of the results. This advanced assessment is introduced in the present study for the subsystem refrigeration machine.

Subsystem refrigeration machine

The optimization potential index for refrigeration machines is [START_REF] Brenner | Optimization potential index (OPI): An evaluation method for performance assessment and optimization potential of chillers in HVAC plants[END_REF]:

OP I RM = 1 - 24 h t=0 h B * el,CP R 24 h t=0 h B el,CP R = 1 - 24 h t=0 h W * el,CP R 24 h t=0 h W el,CP R (2)
where B el,CP R and B * el,CP R is the actual and reference (adequate) electrical exergy input of the compressor, respectively. Since electrical energy is per definition pure exergy, the key figure may be expressed in terms of s u b m i s s i o n p r e p r i n t the actual W el,CP R and reference electrical energy consumption of the compressor W * el,CP R . A quasi-stationary approach is applied for the OPI, where steady-state is assumed over the measurement interval. Consequently, the electrical power consumption is integrated over each measurement interval and then summed up for the daily assessment. For the advanced evaluation, it is proposed to define an additional acceptable boundary given by:

OP I acc RM = 1 - 24 h t=0 h W * el,CP R 24 h t=0 h W acc el,CP R (3)
with W acc el,CP R the acceptable electrical energy consumption of the compressor. As an approach, it is proposed to compute the reference (adequate) and acceptable compressor power consumption with a model of the refrigeration machine (see section 4). Since the present method should be widely applicable in practice, the model should be as general as possible and require few as well as easily obtainable data with state-ofthe-art measuring concepts. Consequently, secondary side temperatures and the evaporator cooling capacity are specified as model input variables, which are commonly measured in field plants (see section 3). The reference compressor electrical power consumption Ẇ * el,CP R is calculated by applying reference secondary side temperatures of the condenser (T * C,in and T * C,out ) together with the cooling load QE and cold water temperatures (T E,in and T E,out ) of the actual situation as model input parameters. In this paper it is proposed to define the reference condenser inlet temperature T * C,in according to:

T * C,in = T C,in -   T C,out -T C,in ln T C,out T C,in -T amb + ∆T HE + T C,out -T C,in 2   (4)
with T C,in and T C,out the secondary side condenser inlet and outlet temperature, respectively. T amb denotes the ambient air temperature and ∆T HE the temperature difference in the cooler heat exchanger, i.e. cooling water outlet to ambient air inlet, according to technical standards. For the latter, a stricter value of 6 K is applied as reported in the technical standard VDMA 24247-8 [START_REF]Energieeffizienz von Klimakälteanlagen. Teil 8: Komponenten -Wärmeübertrager[END_REF] (see Table 1). Similarly, the reference condenser outlet temperature T * C,out is proposed in the present work by:

T * C,out = T C,out -   T C,out -T C,in ln T C,out T C,in -T amb + ∆T HE + T C,out -T C,in 2   (5)
Therefore, these reference temperatures represent adequate values which should be achieved if the hot side hydraulic circuit including the cooler is correctly operated and maintained according to the technical standards.

The acceptable electrical power consumption of the compressor Ẇ acc el,CP R is determined analogously with the refrigeration machine model by applying acceptable secondary side temperatures of the condenser (T acc C,in and T acc C,out ). These temperatures are calculated similarly to Eq. 4 and 5 by applying the less stricter value for ∆T HE according to Table 1.

s u b m i s s i o n p r e p r i n t is integrated to the system and the refrigeration machines as well as the distribution networks (see Fig. 2, subsystem CST) are located underground. The hydraulic circuit supplies seven different buildings with cold water, where the cooling locations (see Fig. 2, subsystem CL) represent air-handling units of ventilation systems in the different buildings with office space cooling as main application. Additionally, three rooftop coolers (see Fig. 2, subsystem DC), 12 circulating pumps and two cold water storages are present in the system. s u b m i s s i o n p r e p r i n t

Investigated modeling approaches

Equation-fit based model

The Comstock model has been chosen as equation-fit based model, as it uses commonly measured quantities in field plants as input variables to determine the compressor electrical power consumption. It is a secondorder polynomial fitting model with 7 fitting parameters a i to determine the compressor electrical power consumption ẆCP R , which has the following form [START_REF] Wang | Empirical model for evaluating power consumption of centrifugal chillers[END_REF]:

Ẇel,CP R = a 1 + a 2 T E,in + a 3 T C,in + a 4 QE + a 5 T E,in QE + a 6 T C,in QE + a 7 Q2 E ( 6
)
with T E,in the evaporator secondary side inlet temperature, T C,in the condenser secondary side inlet temperature and QE the evaporator heat flow rate as input parameters. The coefficients a i reveal no physical characteristic and are determined from measurement data with the lsqcurvefit curve-fitting algorithm preimplemented in MATLAB [10].

Pyhsical lumped parameter model

In the present work, the modified Gordon-Ng model by Folicao et al. [START_REF] Foliaco | Improving the Gordon-Ng Model and Analyzing Thermodynamic Parameters to Evaluate Performance in a Water-Cooled Centrifugal Chiller[END_REF] is applied as physical lumped parameter model, which uses commonly measured quantities in field plants as input variables like the Comstock model. It is a three-parameter model which has the following form [START_REF] Foliaco | Improving the Gordon-Ng Model and Analyzing Thermodynamic Parameters to Evaluate Performance in a Water-Cooled Centrifugal Chiller[END_REF]:

Ẇel,CP R - T C,in -T E,out T E,out QE = a 1 T C,in + a 2 T C,in -T E,out T E,out + a 3 Q2 E + QE Ẇel,CP R T E,out (7)
with T E,out the evaporator secondary side outlet temperature, T C,in the condenser secondary side inlet temperature and QE the evaporator heat flow rate as input parameters. The fitting parameters are determined with the regression routine regress pre-implemented in MATLAB [10], where a 1 denotes the entropy generation in the refrigeration cycle, a 2 the heat losses or gains from the refrigeration machine and a 3 the total heat exchanger thermal resistance.

Refrigeration cycle based model

A simplified refrigeration cycle based model is applied in the present work, which considers the four main refrigeration machine components (evaporator, compressor, condenser and expansion valve) in terms of thermodynamic and heat transfer relations. Fig. 3 shows a schematic of the considered refrigeration cycle as well as the corresponding log(p)-h-diagram with the different refrigerant states. The thermo-physical properties of the refrigerant are computed with the tool REFPROP (Reference Fluid Thermodynamic and Transport Properties Database) [START_REF] Bell | NIST Reference Fluid Thermodynamic and Transport Properties Database (REFPROP) Version 9 -SRD 23[END_REF]. The model uses secondary side temperatures as well as the evaporator cooling load as input variables and incorporates 6 physical parameters (see The evaporator and condenser are modeled with the NTU-ε effectiveness method and are considered as heat exchangers with phase change on the primary side. This represents a special case, where the heat capacity rates of the condensing vapor or evaporating liquid tends towards infinity (nearly isothermal processes). The evaporator effectiveness ε E is then given by [START_REF] Incropera | Principles of heat and mass transfer[END_REF][START_REF] Jin | A parameter estimation based model of water-to-water heat pumps for use in energy calculations programs[END_REF]:

ε E = 1 -e -N T U E ( 8
)
where N T U E denotes the dimensionless parameter number of transfer units of the evaporator, which is defined as [START_REF] Incropera | Principles of heat and mass transfer[END_REF]:

N T U E = (U A) E c p,E ṁE (9)
with (U A) E the overall heat transfer coefficient times the heat transfer area of the evaporator as well as c p,E and ṁE the specific heat capacity and the mass flow rate of the secondary side chilled water, respectively. s u b m i s s i o n p r e p r i n t As the heat capacity and mass flow rate of the secondary side are unknown, their product can be determined with an energy balance over the evaporator. The evaporation temperature of the refrigerant is then given with [START_REF] Incropera | Principles of heat and mass transfer[END_REF][START_REF] Jin | A parameter estimation based model of water-to-water heat pumps for use in energy calculations programs[END_REF]:

T e = T E,in - ∆T E ε E ( 10
)
where T E,in denotes the evaporator secondary side inlet temperature and ∆T E the secondary side inlet to outlet temperature difference. Analogously, Eq. 8 to 10 apply also for the condenser by utilizing the corresponding quantities. With the chosen NTU-ε effectiveness approach, the superheating and subcooling effects in the heat exchangers are not explicitely treated. However, the error is assumed negligible, as it is presumably compensated with the U A parameter when calibrating the model. With the evaporation and condensation temperature, the low and high pressure level, p e and p c (see Fig. 3b), can be determined together with the refrigerant temperature T 1 and T 3 after the evaporator and condenser, respectively. Subsequently, the specific enthalpies and entropies (h 1 , h 3 , s 1 , s 3 ) of the refrigerant can be identified at state 1 and 3. By assuming an isenthalpic expansion process, the specific enthalpy h 4 at state 4 after the expansion valve (see Fig. 3) is equivalent to the specific enthalpy h 3 at state 3. With an energy balance, the evaporator heat flow rate QE is

given by [START_REF] Dinçer | Refrigeration Systems and Applications[END_REF][START_REF] Moran | Fundamentals of engineering thermodynamics: SI version[END_REF]:

QE = ṁr (h 1 -h 4 ) (11)
from which the refrigerant mass flow rate ṁr can be determined. The actual compressor input ẆCP R,act is defined according to [START_REF] Dinçer | Refrigeration Systems and Applications[END_REF][START_REF] Moran | Fundamentals of engineering thermodynamics: SI version[END_REF]:

ẆCP R,act = ṁr (h 2 -h 1 ) (12)
s u b m i s s i o n p r e p r i n t

where the enthalpy h 2 at state 2 after the compressor must be known (see Fig. 3). The compressor isentropic efficiency η isen is given by [START_REF] Dinçer | Refrigeration Systems and Applications[END_REF][START_REF] Moran | Fundamentals of engineering thermodynamics: SI version[END_REF]:

η isen = ẆCP R,isen ẆCP R,act = ṁr (h 2s -h 1 ) ṁr (h 2 -h 1 ) = h 2s -h 1 h 2 -h 1 (13)
which relates the isentropic compressor input ẆCP R,isen with the actual compressor power ẆCP R,act , where h 2s represents the specific enthalpy at state 2s resulting from an isentropic (i.e. reversibel) compression (see Fig. 3). By rearranging Eq. 13, h 2 can be identified, and finally, the compressor electrical power Ẇel,CP R determined with [START_REF] Ozgur | Exergy analysis of refrigeration systems using an alternative refrigerant (hfo-1234yf) to R-134a[END_REF]:

Ẇel,CP R = ẆCP R,act η el,mech (14)
where η el,mech represents the compressor electro-mechanical efficiency to account for mechanical as well as electrical losses. The condenser heat flow rate QC is given with an overall energy balance over the refrigeration cycle [START_REF] Dinçer | Refrigeration Systems and Applications[END_REF]:

QC = QE + ẆCP R,act (15)
The condenser heat flow rate is unknown in the beginning, guessed for the initial iteration and then iteratively computed with the described procedure, until the relative error is lower than 0.1%. A cost function, chosen to be the root-mean-squared error (RMSE, see subsection 4.5) between the modeled and measured compressor electrical power, is minimized to identify the models parameters (see Table 2). These are initially determined with a batch gradient-descent routine [START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF] and then iteratively optimized with the pre-implemented minimization algorithm fmincon in MATLAB [10]. All parameters reveal a physical characteristic and defining them as constants is physically not correct (except in steady-state). However, to the best of the authors knowledge, this simplification is reasonable with the goal of having a widely applicable model which can be employed independently of the compressor and heat exchanger design, even if no technical details of the refrigeration machine are available.

Artifical neural network model

In the present work, a feed-forward ANN model is applied with one input layer, one ouput layer and one hidden layer (see Fig. 4), which has been determined by trial and error. Like the refrigeration cycle based model, the ANN has 5 input parameters (secondary side temperatures as well as evaporator cooling load) and one output parameter, the compressor electrical power Ẇel,CP R . A tan-sigmoid and a linear activation function is used for the hidden and output layer, respectively. The hidden layer consists of 25 neurons (with bias). The amount of neurons was identified by iteratively increasing the number until the performance function, chosen to be the mean-squared error (MSE), converged. The latter is defined as [START_REF] Datta | An optimized ANN for the performance prediction of an automotive air conditioning system[END_REF]:

M SE = 1 n n i=1 ( y i -y i ) 2 ( 16
)
where y i represents the measured value, y i the predicted value and n the number of data points. The number of neurons is substantially lower than the number of training data points, which reduces the risk of overfitting.

s u b m i s s i o n p r e p r i n t The ANN is designed, set-up and trained with the Deep Learning Toolbox (functions feedforwardnet and train)

pre-implemented in MATLAB [10]. One of the most popular learning method is the back-propagation, which is a gradient-descent based procedure [START_REF] Kalogirou | Applications of artificial neural-networks for energy systems[END_REF]. In the present work, the Levenberg-Marquardt back-propagation algorithm is applied for training [START_REF] Hagan | Training Feedforward Networks with the Marquardt Algorithm[END_REF], which is particularly useful for moderate-sized ANNs with the MSE as performance function [START_REF] Hagan | Neural Network Design[END_REF].

Model performance metrics

Many different performance indicators exist, which can be applied to evaluate the model accuracy. In the present work, the following indicators are applied, where y i represents the measured value, y i the predicted value and n the number of data points. One of the most applied index is the root-mean-squared error (RMSE) defined as [START_REF] ¸ahin | Performance analysis of single-stage refrigeration system with internal heat exchanger using neural network and neuro-fuzzy[END_REF]:

RM SE = 1 n n i=1 ( y i -y i ) 2 (17)
The second index is the mean-absolute error (MAE) which has an increased interpretability compared to the RMSE, which is given by [START_REF] Wang | Empirical model for evaluating power consumption of centrifugal chillers[END_REF]:

M AE = 1 n n i=1 | y i -y i | (18)
The coefficient of variation of the root-mean-squared error (CV or CV-RMSE), also denoted relative rootmean-squared error (R-RMSE), indicates if the model has a satisfactory prediction ability, where a small value indicates a high predictive accuracy. The indicator is defined with [START_REF] Wang | Empirical model for evaluating power consumption of centrifugal chillers[END_REF]:

CV = RM SE y i • 100% (19)
s u b m i s s i o n p r e p r i n t where y i is the average of the measured values. Additionally, the coefficient of determination (R 2 ) is applied, which is given by [START_REF] ¸ahin | Performance analysis of single-stage refrigeration system with internal heat exchanger using neural network and neuro-fuzzy[END_REF]:

R 2 = 1 - n i=1 (y i -y i ) 2 n i=1 (y i -y i ) 2 (20)
This indicator relates the sum of the squared residuals (deviation from the predicted and measured value)

with the total sum of squares. The closer the value of R 2 is to 1, the more accurate are the modeled values.

Results and discussion

Refrigeration machine model comparison

Fig. 5 depicts the modeled compressor electrical power consumption of each modeling approach (y-axis) in function of the measured compressor electrical power consumption (x-axis) with a ±10% error band. Table 3 lists s u b m i s s i o n p r e p r i n t which is adequate for a practical application of the model [START_REF] Wang | Empirical model for evaluating power consumption of centrifugal chillers[END_REF][START_REF] Hydeman | Development and Testing of a Reformulated Regression-Based Electric Chiller Model[END_REF]. Consequently, the ANN model is applied for the optimization potential index. By using reference and acceptable temperatures on the condenser secondary side according to technical standards, the reference and acceptable compressor power consumption can be simulated (see subsection 2.2).

Refrigeration machine performance analysis

By applying the method described in subsection 2.2, together with the acquired measurement data from the field plant and the ANN model, the corresponding OPI of each refrigeration machine is determined. The analysis should further demonstrate the usage of the evaluation approach and reveal the performance as well as eventual optimization potentials of the refrigeration machines. According to the available experimental data, refrigeration machine operation occurred mainly from begin of April to end of October. Fig. 6 shows the daily optimization potential index (OPI, y-axis) of the subsystem RM in the field plant under investigation in function of the date (x-axis). The daily OPI is indicated with data points, where the 14-days moving average is represented by a solid line to evaluate the tendency over time. The green, yellow and red zone depicts the adequate, acceptable and inadequate operation condition, respectively. The adequate boundary is always at OPI = 0 according to the key figure definition (comparison of the reference with the actual effort), while the acceptable limit may fluctuate due to dependencies of various parameters at the different operating points (comparison of the reference with the acceptable effort). Additionally, missing data points for the acceptable limit were determined by interpolation, since not all refrigeration machines were running every day.

The optimization potential index of the five different refrigeration machines lies in the range of -0.47 to 0.44 s u b m i s s i o n p r e p r i n t machine 5 yields the highest average key figure among all refrigeration machines of 0.02 (see Fig. 6e). In 44 and 29% of the time in the investigated period the technical requirements are exceeded and fulfilled, respectively.

The refrigeration machine was commissioned on April 13 th , and therefore, no key figures are present before that day. The late commissioning of the machine may also be the reason for the differentiated behavior compared to the other devices. Overall, OP I RM 1 to OP I RM 5 are at least 73% of the time in adequate and acceptable operation range. This indicates a reasonable performance of all refrigeration machines and their hydraulic integration, where mostly little to no optimization potential compared to the state of the art in technology is present.

Interestingly, by examining the moving average OP I avg RM of all refrigeration machines, they reveal a similar operation. Analogue tendencies of an increasing and decreasing OPI are observed. Presumably, this is due to the same type and size of all installed refrigeration machines, e.g. redundancy purposes, where each of them are operated comparably. Moreover, it is revealed that during the warmer months over the year, i.e. June to end of August, all refrigeration machines exceed and fulfill the technical requirements. This leads to the assumption that the chillers are working near to or at the design point and that the hot side hydraulic circuit is properly operating. The refrigeration machines reveal an increase of the OPI in the transition period, i.e.

April, May and September, where mostly an acceptable operation is present. A further noticeable increase in OPI is observed in the colder months October to December, where generally an inadequate operation is present. However, the systems are operating infrequent in this time period and single outliers should not be overly considered. Most likely, the refrigeration machines are not operating at the design point in the mentioned time period or the temperature level on the hot side hydraulic circuit is not ideal, where the latter can result in an increased electrical power consumption of the compressors. According to literature, lowering the condensing temperature (which is highly influenced by the hot side hydraulic circuit) by 1 K, can reduce the energy consumption of the refrigeration plant by up to 2.5% as a rule of thumb [START_REF] Brunner | Klimakälte heute -Kluge Lösungen für ein angenehmes Raumklima[END_REF]. The increased OPI values in the colder months can also be an indicator to make use of free cooling, which is typically active in this period. Furthermore, errors resulting from the ANN model are possible, since errors in part load conditions may have an increased effect. However, to the best of the authors knowledge, these errors are assumed negligible according to the model performance (see subsection 5.1). To finally determine the issues in the subsystem RM in the mentioned time period a detailed analysis would be necessary, assuming the needed experimental data is available.

Conclusions and outlook

The present work demonstrates successfully the use of refrigeration machine models for the integration in the optimization potential index. The practical application is shown by applying the method to a real field installation as a case study. With the introduced additional limits, a more detailed assessment with the evaluation method is achieved and the interpretability of the results is increased. The distinction between adequate (technical requirements exceeded), acceptable (technical requirements fulfilled) and inadequate (potential for s u b m i s s i o n p r e p r i n t improvement) operation according to the state of the art in technology is straightforward and the results are simple to interpret, which is important in practice. In that way, refrigeration plant operators can track the refrigeration machine performance on a daily basis and a simple colored indicator could be realized for the implementation in monitoring systems. To the best of the authors knowledge, an analysis with the described method (revealing possible optimization potentials in a first step), together with a detailed analysis (identify the malfunction in detail in a second step and evaluate if adjustments are worthwhile also from an economic point of view), delivers a target-oriented procedure to analyze the refrigeration plant behavior and to optimize the system efficiency.

Moreover, all examined models reveal an acceptable performance, where the mean-absolute error ranges from 0.8 to 3.4 kW. This values are considered reasonable, with the given range of the measured compressor electrical power consumption. However, only the ANN model reaches coefficient of variation values lower than 5%, which represents an adequate value for a practical application. Consequently, this modeling approach is applied for the determination of the optimization potential index. The investigated refrigeration machines of the field plant generally fulfill or exceed the technical requirements. Refrigeration machine 2 performs best with an average optimization potential index of -0.03. Refrigeration machine 5 performs worst, while being 73% of the time in adequate or acceptable operation. Moreover, all refrigeration machines show potential for improvement in the colder months. This can be an indicator to make use of free cooling, which is typically active in this time period.

As an outlook, future work could cover the investigation of models for the other subsystems in the refrigeration plant. This can increase the level of detail of the analysis and may helps to assess refrigeration plants with a reduced amount of installed measuring equipment. In this context, further measurement data should be collected from other real refrigeration plants, preferably with different cooling capacities and sizes. The practical application of the evaluation method could then be further elaborated in detail and interrelationships between various field plants might be identified.

Acknowledgements

This study was partially funded by the Swiss Federal Office of Energy (SFOE).

s u b m i s s i o n p r e p r i n t

  optimization potential index Mots-clés : méthode d'évaluation du potentiel d'optimisation, modélisation, installation de réfrigération, machine de réfrigération, optimization potential index s u b m i s s i o n p r e p r i n t Highlights Four different refrigeration machine models are evaluated.

Figure 1 :

 1 Figure 1: Optimization potential index scale to determine the operation condition and improvement possibilities with respect to the technical baseline with (a) a basic and (b) an advanced assessment.

Figure 2 :

 2 Figure 2: Simplified piping & instrumentation diagram of a typical refrigeration plant with cold water distribution and a free cooling module, i.e. heat exchanger. Commonly measured and monitored variables of refrigeration machines in field plants are shown in italic.

Figure 3 :

 3 Figure 3: Simplified schematic of (a) the considered refrigeration cycle and (b) the corresponding log(p)-h-diagram with the different refrigerant states 1 to 4 and 2s.

Figure 4 :

 4 Figure 4: Topology of the applied feed-forward neural network with the corresponding input and output variables.

  the values of the different performance indicators for each model with respect to the training / validation and testing data. For the former, measurement data of refrigeration machine 1 (RM1, 12430 data points) installed in the investigated field plant are applied, where a share of 60 and 40% is used for training and validation, respectively. Measurement data of refrigeration machine 2 (RM2, 12250 data points) is applied as testing data set. The measured compressor electrical power ranges from 1 to 195 kW and from 1 to 184 kW in the training / validation and testing data set, respectively, where a large range of different operating conditions are present (part load as well as full load).

Figure 5 :

 5 Figure 5: Comparison of the measured and modeled compressor electrical power consumption of the refrigeration machines in the field plant: (a) equation-fit based model, (b) physical lumped parameter model, (c) refrigeration cycle based model and (d) artifical neural network model.

(

  see Fig.6ato 6e), whereas the maximum difference in OPI of 0.43 between two chillers (RM4 and RM5) is reached on June 15 th . Refrigeration machine 2 yields the lowest average optimization potential index OP I RM 2 of -0.03 (see Fig.6b), while the key figure scatters the least among all refrigeration machines. The system operates approximately 80% of the time below the acceptable boundary OP I acc RM 2 , where 60 and 20% of the time an adequate and acceptable operation of the subsystem is achieved, respectively. Conversely, refrigeration s u b m i s s i o n p r e p r i n t

Figure 6 :

 6 Figure 6: Optimization potential index OP I RM of (a) refrigeration machine 1 to (e) refrigeration machine 5 with adequate (green), acceptable (yellow) and inadequate (red) operation range. The data points (black crosses) represent the daily OPI values and the black solid line indicates the 14-days moving average of the OPI.

Table 1 :

 1 Temperature differences in dry cooler heat exchangers according to VDMA 24247-8[START_REF]Energieeffizienz von Klimakälteanlagen. Teil 8: Komponenten -Wärmeübertrager[END_REF].Fig.2shows schematically the typical refrigeration plant structure with cold water distribution and free cooling together with measured and monitored thermodynamic quantities of refrigeration machines in field plants. As a case study, an existing refrigeration plant installed in the city of Winterthur, Switzerland is investigated in the present work. The field plant includes five refrigeration machines (see Fig.2, subsystem RM) with 950 kW cooling power each and ammonia (R717) as refrigerant. The subsystem refrigeration machine is of particular interest in the present work, where measurement data was collected during the year 2018 with a measuring interval of 5 minutes. Additionally, one free cooling heat exchanger (see Fig.2, subsystem FC)

		Adequate Acceptable
	∆T HE	≤ 6 K	≤ 8 K
	3. Refrigeration plant structure

Table 2

 2

	for details), which are identified from

Table 2 :

 2 Inputs, outputs and parameters of the RC based model.

	Type	Variables
	input	T E,in	evaporator secondary side inlet temperature
		T E,out	evaporator secondary side outlet temperature
		T C,in	condenser secondary side inlet temperature
		T C,out	condenser secondary side outlet temperature
		QE	evaporator heat flow rate
	output	Ẇel,CP R compressor electrical power
	parameter η isen	compressor isentropic efficiency
		η el,mech	compressor electro-mechanical efficiency
		∆T sh	superheating temperature difference
		∆T sc	subcooling temperature difference
		(U A) E	evaporator surface dependent overall heat transfer coefficient
		(U A) C	condenser surface dependent overall heat transfer coefficient

Table 3 :

 3 List of the RMSE, MAE, R 2 and CV values for each modeling approach with respect to the training / validation (RM1) and testing (RM2) data of the refrigeration machines in the field plant.

	Performance indicator	Equation-fit based model	Phyiscal lumped parameter model
		RM1	RM2	RM1	RM2
	root-mean-squared error (RMSE) [kW] 3.23	3.26	3.15	3.20
	mean-absolute error (MAE) [kW]	2.36	2.34	2.44	2.55
	coefficient of determination (R 2 ) [-]	0.995	0.994	0.995	0.995
	coefficient of variation (CV) [%]	6.36	6.34	6.20	6.22
		Refrigeration cycle	Artificial neural
		based model	network model
		RM1	RM2	RM1	RM2
	root-mean-squared error (RMSE) [kW] 4.02	4.24	1.12	1.45
	mean-absolute error (MAE) [kW]	3.16	3.40	0.80	1.07
	coefficient of determination (R 2 ) [-]	0.992	0.991	0.999	0.999
	coefficient of variation (CV) [%]	7.91	8.25	2.21	2.82

Table 3 )

 3 , respectively, where the values are reasonable with the given range of the compressor electrical power. Similar values are achieved with the testing data set, which has not been used for training. The physical lumped parameter (PLP) model reveals a similar performance as the EF model, where the RMSE and MAE values range depending on the data is present in the experimental data of the field plant, which complicates the identification of valid parameters for all load conditions. The predicted power consumption by the artificial neural network (ANN) model shows qualitatively a good agreement with measurements when using both data sets (see Fig.5d). The adequate performance is also demonstrated with R 2 values close to 1 and CV values below 3%, where only approximately 10% of the simulated values are outside of the ±10% error band.When comparing the different models, all approaches perform only slightly worse with the testing data set, where the EF and PLP model reveal a similar performance. This outcome is probably due the large training domain, where all kind of operation conditions are covered. Also, all refrigeration machines are of the same type and size, whereby most likely all of them reveal a similar operating behavior. The RC model underpredicts the compressor power significantly in part load conditions (see Fig.5c), which contributes to the highest key figure values compared to the other investigated models. All performance indicators reveal acceptable values except the coefficient of variation, where only the ANN model reaches values lower than 5%,

	p r e p r i n t
	s u b m i s s i o n

set from 3.15 to 3.2 kW and 2.44 to 2.55 kW, respectively (see Table

3

). Furthermore, the refrigeration cycle (RC) based model performs worst among the investigated modeling approaches. This outcome is most likely due to the constant model parameters, which would vary in reality. A large range of different operating states


Theodore Lytras

Manolis Kogevinas

Hans Kromhout

Anne-Elie Carsin

Josep Maria Antó

Hayat Bentouhami

Joost Weyler

Joachim Heinrich

Dennis Nowak

Isabel Urrutia

Jesús Martínez- Moratalla

José Antonio Gullón

Antonio Pereira Vega

Chantal Raherison Semjen

Isabelle Pin

Pascal Demoly

Bénédicte Leynaert

Simona Villani

Thorarinn Gislason

Cecilie Svanes

Mathias Holm

Bertil Forsberg

Dan Norbäck

Amar J Mehta

Nicole Probst-Hensch

Geza Benke

Kjell Torén

Torben Sigsgaard

Vivi Schlünssen

Mario Olivieri

Paul D Blanc

Roel Vermeulen

Judith Garcia-Aymerich

Deborah Jarvis

Jan-Paul Zock

Occupational exposures and 20-year incidence of COPD: The European Community Respiratory Health Survey

Background: Occupational exposures have been associated with an increased risk of Chronic Obstructive Pulmonary Disease (COPD). However, few studies have related objectively-assessed occupational exposures to prospectively-assessed incidence of COPD, using post-bronchodilator lung function tests. Our objective was to examine the effect of occupational exposures on COPD incidence in the European Community Respiratory Health Survey (ECRHS).

Methods: General population samples aged 20-44 were randomly selected in 1991-1993, and followed up 20 years later (2010-2012). Spirometry was performed at baseline and at follow-up, with incident COPD defined using a Lower Limit of Normal (LLN) criterion for post-bronchodilator FEV1/FVC. Only participants without COPD and without current asthma at baseline were included. Coded job histories during follow-up were linked to a Job-Exposure Matrix, generating occupational exposure estimates to twelve categories of agents. Their association with COPD incidence was examined in log-binomial models fitted in a Bayesian framework.

Findings: 3,343 participants fulfilled the inclusion criteria; 89 of them had COPD at follow-up (1.4 cases/1000 person-years). Participants exposed to biological dust had a higher incidence of COPD compared to those unexposed (RR=1.6, 95% Credible Interval: 1.1 -2.3), as did those exposed to gases & fumes (RR=1.5, 95%CrI: 1.0 -2.2) and pesticides (RR=2.2, 95%CrI: 1.1 -3.8). The combined Population Attributable Fraction for these exposures was 21.0%.

Interpretation: These results substantially strengthen the evidence base for occupational exposures as an important risk factor for COPD.

Key messages

What is the key question?

• What is the effect of occupational exposures on postbronchodilator spirometry-defined COPD incidence?

What is the bottom line?

• Exposure to biological dusts, gases & fumes, and pesticides were associated with increased COPD incidence, and together accounted for 21% of cases in the study population.

Why read on?

• This is the first multicentre prospective study to show an effect of biological dusts and of pesticides on COPD incidence, substantially strengthening the evidence base for occupation as a risk factor for COPD.

Introduction

Chronic Obstructive Pulmonary Disease (COPD) is a progressive respiratory disease characterized by a largely non-reversible obstruction of the airways leading to airflow limitation. [START_REF]Global strategy for the diagnosis, management, and prevention of chronic obstructive pulmonary disease[END_REF] It is a leading cause of morbidity and mortality worldwide, whose impact is expected to further increase as the population ages. [START_REF] Lozano | Global and regional mortality from 235 causes of death for 20 age groups in 1990 and 2010: a systematic analysis for the Global Burden of Disease Study 2010[END_REF][3][START_REF] Herse | Annual costs of chronic obstructive pulmonary disease in Finland during 1996-2006 and a prediction model for 2007-2030[END_REF] Tobacco smoking is the primary risk factor for COPD; [START_REF] Kohansal | The natural history of chronic airflow obstruction revisited: an analysis of the Framingham offspring cohort[END_REF] however, many other environmental factors have been implicated in COPD, including occupational exposures. [START_REF] Eisner | An official American Thoracic Society public policy statement: Novel risk factors and the global burden of chronic obstructive pulmonary disease[END_REF][START_REF] Paulin | Occupational exposures are associated with worse morbidity in patients with chronic obstructive pulmonary disease[END_REF] Numerous studies, both population-based and industry-based, have examined the relationship between occupation and COPD-related outcomes. [START_REF] Eisner | An official American Thoracic Society public policy statement: Novel risk factors and the global burden of chronic obstructive pulmonary disease[END_REF][START_REF] Balmes | American Thoracic Society Statement: Occupational contribution to the burden of airway disease[END_REF][START_REF] Omland | Occupational chronic obstructive pulmonary disease: a systematic literature review[END_REF] An interaction between smoking and occupational exposures has also been observed. [START_REF] Blanc | Occupational exposures and the risk of COPD: dusty trades revisited[END_REF] It has been estimated that about 15% of COPD cases are attributable to exposures at the workplace, [START_REF] Balmes | American Thoracic Society Statement: Occupational contribution to the burden of airway disease[END_REF] with a higher population attributable fraction among non-smokers. [START_REF] Hnizdo | Association between chronic obstructive pulmonary disease and employment by industry and occupation in the US population: a study of data from the Third National Health and Nutrition Examination Survey[END_REF][START_REF] Salvi | Chronic obstructive pulmonary disease in non-smokers[END_REF] Despite a wealth of evidence on COPD risk factors from population-based studies, few such studies have assessed the incidence of COPD in a prospective manner, [START_REF] Terzikhan | Prevalence and incidence of COPD in smokers and non-smokers: the Rotterdam Study[END_REF][START_REF] Afonso | COPD in the general population: prevalence, incidence and survival[END_REF] and even fewer have specifically examined the association between occupational exposures and the incidence of spirometry-defined COPD. [START_REF] Mehta | Occupational exposure to dusts, gases, and fumes and incidence of chronic obstructive pulmonary disease in the Swiss Cohort Study on Air Pollution and Lung and Heart Diseases in Adults[END_REF] The European Community Respiratory Health Survey (ECRHS) is a large multicentre population-based longitudinal study with a long follow-up duration; earlier cross-sectional and longitudinal analyses in this relatively young cohort showed an association between occupation and asthma as well as chronic bronchitis symptoms, but not with accelerated lung function decline or increased COPD incidence. [START_REF] Zock | Occupation, chronic bronchitis, and lung function in young adults. An international study[END_REF]17 The objective of the current analysis was to examine the effect of occupational exposures on post-bronchodilator spirometry-defined COPD incidence in the ECRHS, after 20 years of follow-up.

Methods

ECRHS study overview

The ECRHS is a multicentre longitudinal study initiated in 1991-1993 which enrolled random general population samples aged 20 to 44 years in 55 centres from 23 countries. [START_REF] Burney | The European Community Respiratory Health Survey[END_REF] Participants at baseline (ECRHS I) completed a detailed questionnaire via face-to-face interview and underwent a clinical examination, spirometry and other tests. They were followed-up again between 1998 and 2002 (ECRHS II), and a second time between 2010 and 2012 (ECRHS III). Ethical approval for each centre was obtained from their respective competent bodies, and written informed consent was obtained from all participants.

Study population, spirometry and COPD definition

The population eligible for this study included all participants who completed spirometry at baseline (ECRHS I) and at the second follow-up (ECRHS III). Spirometry was performed according to the ATS/ERS standards for reproducibility, using the maximum value observed per participant for the Forced Volume Capacity (FVC) and Forced Expiratory Volume in 1 second (FEV1). Baseline spirometry at ECRHS I was performed without bronchodilation; spirometry at ECRHS III was performed post-bronchodilation, 15 minutes after administering two 100 mcg puffs of salbutamol using a spacer. We excluded participants with a baseline FEV1/FVC ratio under the Lower Limit of Normal (LLN) for their age, sex and gender according to the GLI-2012 equations. 19 We also excluded participants who reported having current asthma at baseline; current asthma was defined as a positive response to either of the following three questions: "have you had an attack of asthma in the last 12 months?", "are you currently taking any medicines for asthma?" and "have you been woken by an attack of shortness of breath at any time in the last 12 months?". Incident COPD was defined using spirometry only, as a post-bronchodilator FEV1/FVC ratio under the LLN at the end of followup.

Occupational exposure assessment

At both follow-up interviews, participants were asked to provide a detailed list of their occupations and industries from jobs held since the previous study visit that were performed for at least 8 hours a week for at least three months. Each such employment was recorded in free text and subsequently coded in the International Classification of Occupations-88 (ISCO-88) by trained local coders. Occupations were also further grouped into 14 wider job categories, as previously defined. 17 Occupational exposures were assessed by linking the ISCO-88 occupational codes to the semi-quantitative ALOHA(+) Job-Exposure Matrix (JEM). [START_REF] Matheson | Biological dust exposure in the workplace is a risk factor for chronic obstructive pulmonary disease[END_REF] This JEM assigns, for every job code, three grades of exposure (none, low, high) to ten categories of agents (biological dusts, mineral dusts, gases/fumes, herbicides, insecticides, fungicides, aromatic solvents, chlorinated solvents, other solvents, and metals) plus two composites of the above (All pesticides and Vapors/Gases/Dusts/Fumes -VGDF). Exposures for every participant were defined and analyzed for the entire follow-up period, i.e. from ECRHS I to ECRHS III.

Data analysis

Log-binomial regression models were used to estimate the probability of COPD at the end of follow-up period as a function of exposure, covariates and length of follow-up; these models directly provide Relative Risk estimates for each covariate. Given that the incidence of COPD is heavily associated with lung function at baseline, [START_REF] De Marco | Incidence of chronic obstructive pulmonary disease in a cohort of young adults according to the presence of chronic cough and phlegm[END_REF] we further adjusted all models for participants' FEV1/FVC ratio at baseline, expressed as percent predicted according to the GLI-2012 equations. Other covariates used for adjustment were age at baseline, sex, lifetime smoking pack-years, Socioeconomic Status (SES) and early life disadvantage score.

Early life disadvantage score is a composite variable that includes maternal smoking, maternal asthma, paternal asthma, childhood asthma (before age 10), and having a serious respiratory infection before age 5. [START_REF] Svanes | Early life origins of chronic obstructive pulmonary disease[END_REF] SES was defined according to the participants' age of completion of formal education, and classified in three categories: high (>19 years), middle (16-19 years), low (<16 years). We also included quadratic terms for age at baseline and lifetime smoking pack-years, in order to account for any non-linear relationships between these important covariates and COPD incidence. [START_REF] Castaldi | Impact of non-linear smoking effects on the identification of gene-by-smoking interactions in COPD genetics studies[END_REF] Each of the 12 ALOHA(+) exposures was assessed in a univariate fashion, i.e. exposed vs unexposed (to the respective agent). Given the substantial overlap between exposures in certain jobs, we also examined each one in comparison to a common group consisting of those unexposed to all 12 occupational agents under study. In addition for certain exposures we fit bivariate models, i.e. models with two exposures, where the effect of each one is adjusted for the presence of the other. Stratified effects by sex were obtained by including appropriate interaction terms in the models, and dose-response was examined by including separate terms for only low and for ever high exposure. As a sensitivity analysis, we refitted the models after excluding incident asthma cases, i.e. those reporting asthma (as defined above) in either of the two ECRHS follow-up visits. As a further sensitivity analysis, we also fitted "reduced" models without adjustment for SES, early life disadvantage score and baseline FEV1/FVC. A secondary analysis involved using job categories as the exposure, i.e. ever working in a particular category during follow-up, in comparison to a common unexposed group consistently working in white-collar occupations.

All models were fitted in a Bayesian framework with the JAGS software, [START_REF] Plummer | JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling[END_REF] setting non-informative gaussian priors for all fixed-effects parameters, and using 4 chains and 32,000 iterations per chain, discarding the first 2,000 as burn-in. Convergence was checked by visual inspection of the MCMC traceplots and by the Gelman-Rubin statistic. The Bayesian framework avoids convergence issues that have sometimes been reported when fitting log-binomial models in a frequentist setting. To address covariate missingness with respect to lifetime smoking pack-years (for current and ex-smokers), a fully Bayesian imputation sub-model was included, whereby the distribution of pack-years was modelled with a Gamma distribution based on the observed cases, under an ignorable missingness assumption (see Online Supplement).

Based on the Relative Risk estimates, Population Attributable Fractions (PAFs) were calculated both for individual exposures and combinations of exposures (see Online Supplement). [START_REF] Rockhill | Use and misuse of population attributable fractions[END_REF] All analyses were performed with the R statistical environment, version 3.4.1. [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]

Role of the funding source

No sponsor of the study had a role in study design, data collection, data analysis, data interpretation, or writing of the report. The corresponding author had full access to all the data in the study and had final responsibility for the decision to submit for publication.

Results

The flow of ECRHS participants into our study sample is illustrated in Figure 1; in total 3343 participants were analyzed, originating from 24 study centres in 12 countries (Australia, Belgium, Estonia, France, Germany, Iceland, Italy, Norway, Spain, Sweden, Switzerland and United Kingdom). Of these participants, 1409 were never smokers and 1934 were ever smokers (current or ex-smokers at the end of follow-up). The characteristics of our study population are shown on Table 1. Median age at baseline was 34.5 years, and the mean duration of follow-up was 20.0 years (range 18.0-22.4). COPD occured in 96 participants, for an incidence of 1.4 cases/1000 person-years. As expected, COPD incidence in ever smokers was much higher than in never smokers (2.0 vs 0.6 cases/1000 person-years, p<0.001). Lifetime smoking pack-years at followup were missing in 448/1934 ever smokers (23.2%, or 13.4% of the entire study population), for whom Bayesian imputation was performed.

Baseline characteristics were generally similar between the 3343 participants included in the analysis and the 4644 who were eligible (FEV1/FVC > LLN and no current asthma at baseline) but could not be included; those not included in the analysis were slightly more likely to be ever smokers and of lower socioeconomic status than those finally analyzed (Supplementary Table 1).

The percentage of participants exposed ranged across the occupational agents, from a minimum of 4.9% to pesticides to a maximum of 48.4% to VGDF; 1696 participants (50.7%) did not have any low or high occupational exposure during the follow-up period. There was substantial overlap and thus correlation between exposures (Figure 2), particularly between the various types of pesticides and solvents, and also between gases & fumes and other exposures.

The effect (Relative Risk and 95% Credible Interval) of occupational exposures on COPD incidence compared to participants unexposed (to the respective agent) is presented on Table 2. For all exposures except solvents and metals, the crude incidence of COPD among exposed (percentage of participants with COPD at the end of follow-up) was higher than among those unexposed. After adjusting for covariates in log-binomial models, there was a significant effect for biological dust (RR = 1.6, 95% CrI: 1.1 -2.3 vs unexposed) as well as for all pesticides (RR = 2.2, 95% CrI: 1.1 -3.8), and in particular insecticides (RR = 2.3, 95% CrI: 1.1 -4.2); however, the effect for all pesticides and insecticides was based on a small number of cases (10 and 8 respectively). Gases & fumes showed a significant effect compared to participants with no such exposure (RR = 1.5, 95% CrI: 1.0 -2.2), while for mineral dust and VGDF the effect was not significant.

Running the models after excluding incident asthma cases did not materially change the results (Supplementary Table 2); biological dust and gases & fumes showed a marginally stronger effect (RR = 2.0, 95% CrI: 1.2 -3.1, and RR = 1.9, 95% CrI: 1.2 -3.0, respectively), while the effect for pesticides was the same (RR = 1.9, 95% CrI: 0.8 -3.9). We also ran the models using a common, fully unexposed comparator group (Supplementary Table 3); this also did not substantially alter the results, although for gases & fumes the effect fell below the conventional level of statistical significance (RR = 1.4, 95% CrI: 1.0 -2.1).

Including both any biological dust and any pesticide exposure in the same log-binomial model, mutually adjusting each effect for the presence of the other, resulted in a RR of 1.5, 95% CrI: 1.0 -2.2 and 1.8, 95% CrI: 0.9 -3.2 for biological dust and pesticides respectively. Adding also exposure to gases & fumes to this model resulted in increased uncertainty and an RR of 1.3, 95% CrI: 0.8 -2.2 for biological dust, 1.7, 95% CrI: 0.9 -3.2 for pesticides, and 1.2, 95% CrI: 0.7 -2.0 for gases & fumes. Assuming these point Relative Risks and the proportions of exposed participants in our study population, the corresponding Population Attributable Fractions (PAFs) were 10.5% for biological dust, 4.4% for pesticides and 9.2% for gases & fumes; the combined PAF [START_REF] Rockhill | Use and misuse of population attributable fractions[END_REF] for all three exposures was 21.0%.

No differences between men and women were observed in the effect of occupational exposures on the incidence of COPD (Supplementary Table 4). We also ran models for each exposure stratified by intensity, i.e. no exposure vs any low vs any high exposure. In these models there was weak evidence of a dose response relationship for biological dust, but none for pesticides, gases & fumes or other exposures (Table 3).

The additional Relative Risk for high biological dust exposure compared to low exposure was 1.5, 95% CrI: 0.8 -2.8, as there were few incident COPD cases particularly in the high exposure group (10/111, 9.0%).

Among covariates in the models, age and smoking pack-years were strongly associated with COPD, as expected. Men were only slightly more likely than women to suffer from COPD (adjusted RR = 1.2, 95% CrI: 0.8 -1.7). Importantly, there was a strong inverse relationship between %predicted FEV1/FVC at basline and COPD incidence at the end of the follow-up period (Figure 3); each percentage point of increase lowered the incidence of COPD by a factor (RR) of 0.8, 95% CrI: 0.8 -0.9. Early life disadvantage score did not affect COPD incidence (RR per unit change = 1.0, 95% CrI: 0.8 -1.3), and neither did medium (RR = 0.9, 95% CrI: 0.6 -1.3) or low socioeconomic status (RR = 0.8, 95% CrI: 0.4 -1.5).

Associations between occupational exposures and COPD incidence were not substantially modified by omitting adjustments for SES, early life disadvantage score and/or %predicted FEV1/FVC at baseline (Supplementary Table 5).

The secondary analysis by job category indicated a higher incidence of COPD in participants who had worked in the transport industry (1.7, 95% CrI: 0.9 -3.2), in the wood, paper and textile industry (2.1, 95% CrI: 0.9 -4.5) and in agriculture, fishery and forestry (2.3, 95% CrI: 0.9 -5.0); however the precision of these estimates was limited by the low number of study participants working in each job category (Supplementary Table 6). Occupations belonging to these three categories were very frequently exposed to biological dust, pesticides and gases/fumes. Overall the most frequent occupations with exposure to these agents are shown on Supplementary Table 7.

Discussion

Our study found that exposure to biological dust in a general population cohort did result in an increased incidence of spirometrically-defined COPD over the course of 20 years. Biological dust has been shown to affect lung function decline and chronic respiratory symptoms in several worker-based studies such as in farmers, cotton textile and woodworkers. [START_REF] Christiani | Longitudinal changes in pulmonary function and respiratory symptoms in cotton textile workers. A 15-yr follow-up study[END_REF][START_REF] Jacobsen | Cross-shift and longitudinal changes in FEV1 among wood dust exposed workers[END_REF][START_REF] Kahraman | Longitudinal pulmonary functional loss in cotton textile workers: a 5-year follow-up study[END_REF][START_REF] Mandryk | Work-related symptoms and dose-response relationships for personal exposures and pulmonary function among woodworkers[END_REF] An association with COPD prevalence has also been demonstrated in a community-based cross-sectional study, [START_REF] Matheson | Biological dust exposure in the workplace is a risk factor for chronic obstructive pulmonary disease[END_REF] and recently in another one from Australia. [START_REF] Alif | Occupational exposure to pesticides are associated with fixed airflow obstruction in middle-age[END_REF] To our knowledge this is the first study to demonstrate an effect of biological dust exposure on the incidence of COPD in a prospective fashion in a general population cohort. This effect remained strong even after adjustment for other exposures or exclusion of incident asthma cases, and also exhibited some evidence of a dose-response relationship, i.e. a higher incidence of COPD with high exposure to biological dust. As such the study is a significant addition to the evidence base about the role of biological dust on respiratory health.

Exposure to pesticides has also been linked with asthma and possibly COPD in a variety of studies. [START_REF] Alif | Occupational exposure to pesticides are associated with fixed airflow obstruction in middle-age[END_REF][START_REF] Doust | Is pesticide exposure a cause of obstructive airways disease?[END_REF][START_REF] Negatu | Occupational pesticide exposure and respiratory health: a large-scale cross-sectional study in three commercial farming systems in Ethiopia[END_REF] More recently an association of pesticide exposure with accelerated lung function decline was observed in a longitudinal community-based study from the Netherlands; this association was stronger for smokers. [START_REF] De Jong | Association of occupational pesticide exposure with accelerated longitudinal decline in lung function[END_REF] Our study adds further prospective confirmation that pesticides increase the incidence of COPD. Interestingly in our study exposure to pesticides occured in the 1990s and 2000s, while in the Dutch study exposure happed mostly before 1990; this may indicate that contemporary changes in active ingredients, pesticide application methods and personal protective measures may have been insufficient to mitigate the risk to workers.

On the other hand, the lack of association between mineral dust exposure and COPD incidence may be explained by the composition of jobs included in this category (Supplementary Table 7). Older studies found a link between mining (especially coal mining), as well as various industrial occupations, and COPD. Such mining and industrial jobs were a minority in our study population, which was recruited in the early '90s.

Instead, the most frequent jobs with exposure to mineral dust in our study were truck and lorry drivers for men, and cleaners for women.

Disentangling the effects of multiple overlapping occupational exposures is a particular challenge. Our approach was twofold, i.e. we compared each exposure with those unexposed, and separately with a common fully occupationally unexposed comparison group. The first choice of comparator is simpler, but includes other occupational exposures which may attenuate any association with the exposure under study. The second comparator (fully unexposed group) avoids this problem, but may increase the potential for residual confounding particularly with respect to SES, as the people in this group (mostly white-collar workers) are likely to be different in more ways than a model can handle. Nevertheless, in our case we found similar results with both approaches. We then followed up with bivariate and multivariate models, i.e. models that include more than one occupational exposure, as well as an analysis by job category; however this requires a large sample size to maintain precision in the effect estimates.

Occupational exposures have been associated with both asthma and COPD incidence in population studies; [START_REF] Mehta | Occupational exposure to dusts, gases, and fumes and incidence of chronic obstructive pulmonary disease in the Swiss Cohort Study on Air Pollution and Lung and Heart Diseases in Adults[END_REF][START_REF] Kogevinas | Exposure to substances in the workplace and new-onset asthma: an international prospective population-based study (ECRHS-II)[END_REF] COPD and asthma are two different disease entities, modulated by distinct pathophysiological mechanisms, but with considerable clinical and symptom overlap. [START_REF] Postma | The Asthma-COPD Overlap Syndrome[END_REF] As a result, in our basic analyses we did not exclude incident asthma cases. A sensitivity analysis with incident asthma cases excluded did not substantially modify the results.

In all models we adjusted for the FEV1/FVC ratio at baseline, expressed as %predicted, which we found to be a strong predictor of subsequent COPD incidence. This is hardly surprising, since COPD is not a stochastic event, but one that results from a progressive lung function decline below a given threshold. As such the "distance to be covered" between baseline FEV1/FVC and the LLN threshold is critical, and we believe it should be adjusted for in future prospective studies of COPD incidence.

The study has several strengths. Being a population-based study increases generalizability of the findings, especially given the prospective design and follow up of 20 years, one of the longest to date. Complete occupational histories were collected, and a full spectrum of exposures for the follow-up period was determined using a JEM instead of self-report, which could be more vulnerable to bias. We were thus able to explore the effect of particular exposures on COPD incidence and their contribution to the disease burden.

Confounders were tightly controlled: we included adjustments for SES and baseline lung function, and we accounted for non-linear relationships with age and smoking pack-years using quadratic terms; therefore residual confounding by intensity of smoking is most unlikely. Lifetime pack-years of smoking were missing in many current or ex-smokers of our study population; we used fully Bayesian imputation to handle the problem and draw reliable inference while reflecting the appropriate uncertainty, assuming ignorable missingness (see Online supplement).

On the other hand, COPD incidence in our study population was much lower than in other published studies; [START_REF] Terzikhan | Prevalence and incidence of COPD in smokers and non-smokers: the Rotterdam Study[END_REF][START_REF] Pallasaho | Combined effect of smoking and occupational exposure to dusts, gases or fumes on the incidence of COPD[END_REF][START_REF] Lindberg | Ten-year cumulative incidence of COPD and risk factors for incident disease in a symptomatic cohort[END_REF] the most likely reason for this was the still fairly young age of our cohort at follow-up (median age 55, range 39 -68 years). Also the use of pre-bronchodilation spirometry at baseline to filter out prevalent COPD cases, may have led to the exclusion of some participants with reversible or borderline obstruction at baseline. Because of the low number of incident COPD cases, especially in lifetime non-smokers, we could not stratify our analysis by smoking status to check for effect modification. For the same reason several of our effect estimates were characterized by low precision, especially those stratified by sex and by intensity of exposure. In addition, we defined COPD with post-bronchidilation spirometry only and not based on symptoms of chronic bronchitis, which would have raised specificity but severely limited the number of incident cases in this analysis. Another limitation stems from the fact that the outcome of COPD was only assessed at the end of follow-up, thus could in theory have predated the exposure. However, out of 56 occupationally exposed participants with incident COPD, most were already working at an exposed job at the start of follow-up, or had a pre-bronchodilation FEV1/FVC>LLN at ECRHS II after having started working at exposed jobs. Only for 7/56 participants (12.5%) it could not be disproven that COPD might have occured before the occupational exposure. Still this is very unlikely, as COPD incidence is strongly related to age, and because participants developing COPD might be less likely to start working in exposed jobs unless previously exposed.

Nevertheless, the fact that some occupational exposure effects (namely biological dust, gases & fumes and pesticides) were observed even in this relatively young study population, is highly important from a public health point of view. A large proportion of workers have these exposures (49% among ever smokers in our study population), and their association with COPD incidence translates into a substantial number of attributable cases in the population (an estimated 21.0% in this analysis); this suggests that up to one in five new COPD cases among middle-aged people in Western countries could be prevented by avoiding or controlling occupational exposures in contemporary jobs. Our study confirms previous findings regarding the role of occupation on the burden of COPD, [START_REF] Eisner | An official American Thoracic Society public policy statement: Novel risk factors and the global burden of chronic obstructive pulmonary disease[END_REF][START_REF] Balmes | American Thoracic Society Statement: Occupational contribution to the burden of airway disease[END_REF] and provides strong prospective evidence and direct quantitative estimates for the impact of occupational exposures on COPD incidence. Future studies should clarify whether these effects are modified by smoking, their interplay with asthma, and further detail the risks involved with respect to particular occupations, activities and noxious agents. 1.0 (0.5 -1.9) 0.7 (0.2 -1.7) Relative Risks adjusted for sex, age, pack-years of smoking, FEV1/FVC ratio at baseline (%predicted), socioeconomic status and early life disadvantage score.

Figure 1 :Figure 2 :

 12 Figure 1: Flow chart of ECRHS participants into the study population, and reasons for exclusion

Figure 3 :

 3 Figure 3: Percentage of incident COPD cases according to %predicted FEV1/FVC at baseline, observed and model-fitted. For the observed cases, the percentage is calculated considering the participants with %predicted FEV1/FVC ±1 percentage point.

Table 1 : Characteristics of the study population

 1

		Total	Ever smokers	Never smokers	Men	Women
	Number of participants	3343	1934	1409	1638	1705
	Person-years of follow-up	66854	38648	28207	32709	34145
	Mean follow up &	20.0	20.0	20.0	20.0	20.0
	ramge (y)	(18.0 -22.4)	(18.0 -22.4)	(18.2 -22.1)	(18.0 -22.3)	(18.2 -22.4)
	Number of
	incident COPD	96	79	17	51	45
	cases
	COPD incidence
	per 1000 person-	1.4	2.0	0.6	1.6	1.3
	years
	% men	49.0	51.7	45.4	-	-
	Mean age at baseline (y)	34.2	34.7	33.5	34.3	34.0
	Mean lifetime
	pack-years of smoking at	10.6	20.6	-	13.3	7.9
	follow-up
	% current asthma at follow-up	8.1	8.4	7.7	6.2	10.0
	% of participants exposed (low or high)
	Biological dust	32.3	32.5	32.0	30.0	34.4
	Mineral Dust	25.2	28.1	21.2	36.7	14.1
	Gases & fumes	43.5	45.8	40.4	51.7	35.7
	Vapors, Gases, Dusts & Fumes	48.4	50.1	46.1	56.4	40.7
	Herbicides	2.2	2.1	2.4	3.2	1.3
	Insecticides	3.4	3.5	3.3	4.9	1.9
	Fungicides	3.9	4.2	3.5	6.3	1.7
	All pesticides	4.9	5.4	4.3	7.7	2.2
	Aromatic solvents	16.4	17.1	15.4	26.7	6.4
	Chlorinated solvents	13.1	14.2	11.7	21.1	5.5
	Other solvents	28.5	27.8	29.4	32.1	25.0
	Metals	12.6	13.6	11.2	22.9	2.7

Table 3 : Associations between occupational exposures and COPD incidence, stratified by intensity of exposure

 3

					Relative Risk	Relative Risk
					(95% Credible	(95% Credible
					Interval),	Interval),
		Cases in	Cases in	Cases in	ever low vs no	ever high vs no
		unexposed	ever low-exposed	ever high-exposed	exposure	exposure
	Biological dust	55/2264 (2.4%)	31/905 (3.4%)	10/174 (5.7%)	1.5 (1.0 -2.2)	2.2 (1.2 -4.0)
	Mineral Dust	65/2501 (2.6%)	18/554 (3.2%)	13/288 (4.5%)	1.0 (0.6 -1.7)	1.3 (0.7 -2.2)
	Gases & fumes	41/1888 (2.2%)	45/1102 (4.1%)	10/353 (2.8%)	1.6 (1.1 -2.5)	1.0 (0.5 -1.9)
	Vapors, Gases,	40/1725	35/1067	21/551	1.3	1.3
	Dusts & Fumes	(2.3%)	(3.3%)	(3.8%)	(0.8 -2.0)	(0.8 -2.3)
	Herbicides	91/3269 (2.8%)	1/49 (2.0%)	4/25 (16.0%)	0.5 (0.0 -3.0)	2.8 (1.0 -5.9)
	Insecticides	88/3229 (2.7%)	3/54 (5.6%)	5/60 (8.3%)	2.0 (0.5 -5.2)	2.2 (0.8 -4.5)
	Fungicides	88/3211 (2.7%)	4/72 (5.6%)	4/60 (6.7%)	1.5 (0.5 -3.7)	2.0 (0.7 -4.6)
	All pesticides	86/3179 (2.7%)	5/94 (5.3%)	5/70 (7.1%)	2.0 (0.7 -4.4)	2.1 (0.8 -4.3)
	Aromatic solvents	80/2796 (2.9%)	14/495 (2.8%)	2/52 (3.8%)	0.8 (0.5 -1.4)	1.3 (0.2 -4.6)
	Chlorinated	83/2904	9/320	4/119	0.8	0.7
	solvents	(2.9%)	(2.8%)	(3.4%)	(0.4 -1.6)	(0.2 -1.7)
	Other solvents	71/2391 (3.0%)	22/883 (2.5%)	3/69 (4.3%)	0.8 (0.5 -1.2)	1.4 (0.4 -4.0)
	Metals	82/2922 (2.8%)	10/287 (3.5%)	4/134 (3.0%)

Funding

Financial support for ECRHS I, for the local centres included in this study:

Australia: Asthma Foundation of Victoria, Allen and Hanbury's, Belgium: Belgian Science Policy Office, National Fund for Scientific Research, Estonia: Estonian Science Foundation, grant no 1088, France: Ministère de la Santé, Glaxo France, Insitut Pneumologique d'Aquitaine, Contrat de Plan Etat-Région Languedoc-Rousillon, CNMATS, CNMRT (90MR/10, 91AF/6), Ministre delegué de la santé, RNSP, France; GSF, Germany: Bundesminister für Forschung und Technologie, Italy: Ministero dell'Università e della Ricerca Scientifica e Tecnologica, CNR, Regione Veneto grant RSF n. 381/05.93, Norway: Norwegian Research Council project no. 101422/310; Spain: Fondo de Investigación Sanitaria (#91/0016-060-05/E, 92/0319 and #93/0393), Hospital General de Albacete, Hospital General Juan Ramón Jiménez, Dirección Regional de Salud Pública (Consejería de Sanidad del Principado de Asturias), CIRIT (1997 SGR 00079) and Servicio Andaluz de Salud; Sweden: The Swedish Medical Research Council, the Swedish Heart Lung Foundation, the Swedish Association against Asthma and Allergy; Switzerland: Swiss national Science Foundation grant 4026-28099; UK: National Asthma Campaign, British Lung Foundation, Department of Health, South Thames Regional Health Authority.

Financial Support for ECRHS III: Australia: National Health & Medical Research Council, Belgium: Antwerp South, Antwerp City: Research Foundation Flanders (FWO), grant code G.0.410.08.N.10 (both sites), Estonia: Tartu-SF0180060s09 from the Estonian Ministry of Education. France: (All) Ministère de la Santé. Programme Hospitalier de Recherche Clinique (PHRC) national 2010. Bordeaux: INSERM U897 Université Bordeaux segalen, Grenoble: Comite Scientifique AGIRadom 2011. Paris: Agence Nationale de la Santé, Région Ile de France, domaine d'intérêt majeur (DIM) Germany : Erfurt: German Research Foundation HE 3294/10-1 Hamburg: German Research Foundation MA 711/6-1, NO 262/7-1 Iceland: Reykjavik, The Landspitali University Hospital Research Fund, University of Iceland Research Fund, ResMed Foundation, California, USA, Orkuveita Reykjavikur (Geothermal plant), Vegagerðin (The Icelandic Road Administration (ICERA). Italy: All Italian centres were funded by the Italian Ministry of Health, Chiesi Farmaceutici SpA, in addition Verona was funded by Cariverona foundation, Education Ministry (MIUR). Norway: Norwegian Research council grant no 214123, Western Norway Regional Health Authorities grant no 911631, Bergen Medical Research Foundation. Spain: Fondo de Investigación Sanitaria (PS09/02457, PS09/00716 09/01511) PS09/02185 PS09/03190), Servicio Andaluz de Salud, Sociedad Española de Neumología y Cirurgía Torácica (SEPAR 1001/2010); Sweden: All centres were funded by The Swedish Heart and Lung Foundation, The Swedish Asthma and Allergy Association, The Swedish Association against Lung and Heart Disease. Fondo de Investigación Sanitaria (PS09/02457 Barcelona: Fondo de Investigación Sanitaria (FIS PS09/00716) Galdakao: Fondo de Investigación Sanitaria (FIS 09/01511) Huelva: Fondo de Investigación Sanitaria (FIS PS09/02185) and Servicio Andaluz de Salud Oviedo: Fondo de Investigación Sanitaria (FIS

Author Contributions