Similarly, birdsfoot trefoil and sainfoin were found to inhibit the proteolytic bacterial population [49]. A study analysing the effect of different tannin sources on CH4 emissions found that CT-rich (acacia and quebracho tannins) and HT-rich (chestnut and valonea tannins) affect rumen fermentation differently. At concentrations above 5% DM, in addition to a significant decrease in CH4 emissions, there was also a negative effect on total VFA production. CT-rich extracts reduced the acetate/propionate ratio significantly at a concentration higher than 5%. However, the ratio was not affected by HT extracts, indicating that they had a stronger impact on methanogen population in comparison with substrate fermentation. Only valonea extracts (5% *w*/*w*) were able to reduce CH4 emissions without any negative impact on fermentation and VFA profile. This indicates that classification based solely on tannin concentration or the type of tannins (HTs and CTs) present in the feed is not sufficient to determine their potential to reduce CH4 emissions [86]. Similarly, a study was performed on CT-rich forages from Texas to determine the effect of different functional features (PPC and antioxidative activity) on CH4 emissions. No effect of PPC was found on CH4 abatement, whereas the correlation between antioxidative property of tannins and decrease in CH4 emissions was significant. In contrast to previous studies, the decrease in acetate/propionate in this study was not correlated with a decrease in CH4 emissions [90]. The results from these studies further reinforce the need for CT structural characterization in addition to concentration, in order to make an accurate assessment of their impact on ruminant nutrition.

### **4. Existing Research Gaps and Future Directions**

As discussed in the previous sections, several studies have tried to explore the properties that affect tannin astringency. However, the variations in the results obtained, and their lack of reproducibility, hinder their field-scale applicability. Furthermore, their structural complexity and the varied forage chemical composition among different species present difficulties for understanding the implications that CTs have for ruminant nutrition and particularly their antimethanogenic potential. Although several studies have identified a large variability in both the concentration and structure of CTs across species and their cultivars [22,32,74], few studies have analysed the implications of this variability on the observed bioactivity. In the following sections, we present a brief overview about the factors responsible for current situation, with an apparent incongruity regarding the influence of tannins on ruminant nutrition. We also discuss the frequently used analytical techniques for qualitative and quantitative analysis of tannins and the underlying problems associated with them. Our aim here is to illustrate the importance of optimized tannin analyses and inclusion of tannin structural features in animal studies to overcome inconsistent animal responses. By avoiding these factors, which cause substantial variation in the reported studies, we can focus more precisely on CT–animal interactions.

### *4.1. Experimental and Analytical Incongruities*

Tannin concentration and composition in plant has been reported to be substantially influenced by changes in environmental conditions, as well as by plant species and its phenological stage [40]. The preparation and handling of tannin extracts can also cause alterations in quantification of tannins [121]. In order to ensure accurate determination of tannin concentration and composition in a plant tissue, handling and storage protocols should be followed, as CT concentration is highly influenced by the environmental factors. Quantitative analyses are essential for determining CT bioavailability in the samples and spectrophotometric assays are routinely utilized due to their rapid and low-cost analysis. Due to their structural complexity, the number of derivatisation and analytical techniques are few and they have certain limitations. Substantial information about the activity of CTs can be obtained by analysing specific structural traits of tannins as it is difficult to isolate individual large polymeric CT units compared to dimers or trimers [122]. The complexity of CT structures means that they are frequently analysed by a method where multiple techniques with different functions are integrated together.

### 4.1.1. Growth Conditions of Experimental Plants

Tannin concentrations in plants can be up to 20% of their total dry weight [38]. Although, CTs are found in different parts of the plant and they are predominantly concentrated in young leaves and flowers [105,123]. The concentration of tannins in tropical plants is, on average, higher than in temperate plants, yet there is substantial variation across seasons and environmental factors [62]. Drought, nutrient availability, and other conditions during plant growth have also been shown to affect CT concentration and composition. Although the effect of these abiotic stresses on the CT composition has not been well researched [27,28,124], they have been shown to produce incremental effects in CT concentrations [84,125]. Accordingly, it was observed that the concentration of CTs in sulla was higher in the summer than in spring [54]. Similarly, *Quercus rubica* had higher tannin concentration and less polymerized tannins when grown in dry conditions compared with wet conditions [62]. Thus, it is important to account for and report the precise experimental conditions because of their potential to affect the observed tannin concentration, composition, and bioactivity [126].

### 4.1.2. Sample Preservation and Storage

After harvest, sample preservation plays an important role in the quantification of CTs, as their extraction and quantification are heavily influenced by biotic factors. For precise tannin concentration and composition analysis, samples should be freeze-dried, rather than air or oven dried. In a direct comparison, hay-drying of samples from purple prairie clover (*Dalea purpurea* Vent.) resulted in a slight decrease in ECTs from 70.2 to 64.1 g kg DM−1 , while the protein-bound CTs increased from 9.0 to 12.4 g kg DM−1 . With ensiled samples, the differences were even more pronounced, and in these samples, ECTs decreased to 27.4 g kg DM−1 , while protein-bound CTs increased their concentration to 44.3 g kg DM−1 [127]. Thus, while the total CT concentration did not differ, without freeze drying at least a part of the CTs, it can change from the available form to the protein-bound form, which is often not accounted for in the studies analysing bioactivity of tannins. Similarly, when comparing different drying methods for the concentration of HTs in white birch (*Betula pubescens*), oven drying reduced the ET concentration significantly from 10.9 to 8.4 g kg DM−1 , while simultaneously increasing the concentration of insoluble ET from 0.8 to 2.4 g kg DM−1 . In this study, neither storing the samples at −20 ◦C for 3 months prior to drying nor vacuum or air drying resulted in a decrease in ellagitannins or total HTs despite a minor but non-significant decrease in air dried samples [121]. This indicates that the adverse effect from air and oven drying increases with the air temperature, yet the short-term effect at room temperature appears to be negligible. The effect of temperature during post drying storage is less clear, and storage at 25 ◦C for three weeks reduced tannins in walnuts by 20–40% (dependent on the subsequent extraction technique), with large parts of the reduction having occurred in the first week [128]. However, with regard to the HTs from birch, a one-year storage period at room temperature (22 ◦C) yielded lower tannin concentrations compared to samples that were stored in a freezer, with 17.4 and 19.1 g kg DM−1 , respectively. This value was still higher than the HT concentration in the samples stored in a refrigerator at 4 ◦C, which yielded 15.7 g kg DM−1 [121]. Contrasting results were obtained by Kardel [129], where samples stored in a refrigerator for one year had on average 3.5% higher CT concentrations than samples stored at room temperature. In this study, also storing the samples in an oven at 60 ◦C for 5 days had no measurable impact on CT concentrations.

### 4.1.3. Tannin Extraction

While the vast majority of studies use aqueous acetone to extract tannins of any kind, with either 70 or 80% acetone, there is no clear indication about the superiority of any extraction method yet. Some studies have also used methanol or hot water extraction as well. Accordingly, the extraction yields with hot water were 6% higher on average, compared to a water/methanol (1:1) extraction solvent, and 13% higher compared to

an acetone/water/formic acid (70:29.5:0.5) solvent [129]. This study only determined the tannin concentration and did not evaluate potential changes in the structure due to the high temperatures. Contrary to these findings, Salminen [121] found that aqueous acetone extraction yielded on average 41% more extractable HT, compared to aqueous methanol extraction, although the study did not test hot water extraction. Pure acetone, however, yielded the lowest HT concentrations, with a reduction of almost 75% compared to acetone/water (70:30). The extraction yield of the 70% acetone was increased even further by 29%, if ascorbic acid was added to the acetone water mixture, presumably because it prevented oxidation of the HTs [121]. This is in accordance with the findings of Chavan [130] and Hagerman [131] for CTs, where 70% acetone also provided the highest extraction yields compared to all methanol mixtures and acetone mixtures with higher acetone concentrations. Acidification of aqueous acetone extraction solvent with 1mL of concentrated HCl further increased the extraction yields by around 10%.

Some recent studies have, however, indicated in general much lower performance from maceration-based techniques compared to techniques such as ultrasonic baths, microwave-assisted extraction (MAE), and Soxhlet extraction. Aspé [132] identified much larger cell wall destruction from these last three techniques, which resulted in generally much higher extraction yields compared to maceration techniques, where only minor cell wall damage has occurred. According to Chupin, et al. [133], MAE effectiveness depends on the particle size and increases with small particles. However, they generally did not identify an effect of MAE on the structure of the tannins.

In addition to the solvent, the extraction conditions can also affect the efficiency of extraction. Extraction at 4 ◦C in darkness led to around 14–17% higher recovery of CTs compared to extraction at room temperature, also in dark conditions [134].

### 4.1.4. Major Analytical Techniques I: Spectrophotometric Assays

Spectrophotometric assays are routinely used in studies due to their rapid and lowcost analysis. Vanillin and HCl-butanol assays are conventionally employed to quantify CTs. Both these assays have been widely used in the majority of studies due to their accessibility and low cost for quantifying CTs. In principle, results from the HCl-butanol assay are more reproducible than the vanillin assay. Vanillin assay with methanol has been found to lack specificity as in addition to CTs, it complexes with flavan-3-ols and dihydrochalcones [84,135]. In the HCl-butanol assay, CTs are depolymerized oxidatively to form bright coloured anthocyanidins in the presence of mineral acids [135]. The HCl- butanol assay in some instances lacks specificity and leads to over- or underestimation of results as the colour is dependent on the interflavonoid linkages and 5-OH groups [126]. For accurate quantification by the HCl-butanol assay, freeze-dried samples should be preferred as under heat treatment, tannins can bind to other macromolecules which could undergo oxidation. Macromolecules containing phenolic groups can be oxidised to form quinones which can lead to condensation reactions between tannins and other macromolecules, thereby preventing the release of anthocyanidins from the modified tannins [84]. Furthermore, the addition of iron increases the specificity of this assay [135]. The addition of acetone to the assay also leads to the complete dissociation of CTs from the plant material and this inclusion has further refined this method [136]. A wide range of studies use tannic acid, catechin, or leucocyanidin as standards for these assays. This is also one of the sources for which there is a huge variation in the results, as seen in Table 2 [84,129].

4.1.5. Major Analytical Techniques II: Liquid Chromatography Coupled with Mass Spectrometry

High-Performance Liquid Chromatography (HPLC) has been proven to be a competent and rapid method for the analysis of polyphenols apart from the highly polymerized oligomers [48,150]. Reverse phase LC (RPLC) is a commonly used chromatographical technique to analyse CTs, ranging from monomers to tetramers and in some cases, their isomers distinctly [151]. With increased CT polymerization, the intelligibility of the chromatogram in RPLC decreases due to the presence of unresolved peaks. The combination of fluorescence detection with RPLC leads to increases in selectivity and sensitivity of the method [152]. A UV-DAD detector is most frequently complemented with LC for the determination of CTs. It also helps in the direct classification of polyphenols into different subgroups such as flavonoids, ellagitannins, gallic acid derivatives, and caffeic acid derivatives, etc. [153]. Recently, to increase the specificity and resolution of the LC analysis, separation techniques have been coupled with ESI-MS or matrix-assisted laser desorption/ionization time-of-flight (MALDI-TOF) MS [154]. These are used extensively for the analysis of CTs in plant and food material such as cocoa, grapes, wine, and birch species, etc. [155,156]. Soft ionization methods such as ESI or MALDI are used to ionize non-volatile analytes such as biopolymers and detect highly polymerized CTs. MALDI-MS has identified procyanidins of degree of polymerization (DP) of 15 in unripe apples, for seed coats of soyabeans until DP of 30 [157]. Similarly, the combination of Hydrophilic Interaction Chromatography (HILIC) x RPLC with fluorescence detection and electrospray full scan mass spectrometry (ESI-MS) resulted in high resolution analysis. This method was able to detect procyanidins with DP value of 16 and gallolylation degree of 8 in the grape seed extracts [158]. These methods are constantly evolving and are now able to provide rapid quantitative and qualitative results. One such method is the Engstrom method which utilizes ultra-high performance liquid chromatography (UPLC) separation coupled with DAD and negative ion ESI-MS to generate a polyphenolic profile directly from plant extracts. In addition to the quantification of different polyphenolic groups, it provides an insight into the composition of flavonols, CTs, HTs, and structural features of condensed tannins [159,160]. These methods provide a great deal of information on tannin structural diversity, but due to high operational costs, their use is not yet widespread.

### *4.2. Influence of CT Structural Features on Ruminant Nutrition Is Still Ambiguous*

The additive influence of CT structural features on the antimethanogenic potential of CT forages remains largely unexplored. Although research on CT structure and functional features has progressed immensely [64,158,160], only a small number of studies have assessed their impact on CH4 abatement. These studies have shown that CT composition, and the concentration of CTs present in forages, are both significant determinants of their

antimethanogenic potential [24,90,95,161]. One of the major reasons is that CTs exist as highly polymerized structures, so large polymers cannot be easily purified as individual compounds and are studied in terms of certain structural features such as molecular weight, polymer size, and prodelphinidin proportion [33]. High structural variability across and within the species adds to the difficulty for assessing their structure–activity relationship. This has been shown in studies on sainfoin cultivars where antimethanogenic potential was found to be highly variable. Hatew [24] studied the intraspecies variability by analysing 46 different accessions of sainfoin with CT concentrations ranging from 0.6 to 2.8% of DM for their antimethanogenic potential. Emissions were analysed based on CT structural properties, i.e., mDP (12 to 84), percentage of trans isomers (12 to 34%), and PD (52.7 to 94.8%) in CT. These properties have been associated with the astringency of CTs. It was observed that PD percentage was a primary CT structural characteristic responsible for reducing CH4 emissions in this in vitro study [24]. Weight-average Mw of CTs was shown to have little impact on reduction in CH4 production (R2 = 0.0009) from North American native forage plants [117]. Nevertheless, it is important to note that the variation in CH4 reduction potential of different species also arises from plant morphology, CT interaction with feed components, and the presence of other plant secondary metabolites [162,163]. The impact of CT structural features was found to be more pronounced in the studies conducted on CT extracts from plants. The additive effect of other secondary metabolites and forage quality parameters could be voided by the addition of purified CT extracts in the feed. Studies have shown that inclusion of CT extracts (40 mg/g DM) from leucaena (hybrid-Bahru) and mangosteen (*Garcinia mangostana* L) peel could reduce CH4 emissions by 45 and 35 percent, respectively. In both studies, *Panicum maximum* substrate was used as control. The inhibitory effect of mangosteen peel extracts (Mw = 2081) was milder than leucaena (Mw = 2737) owing to its lower MW but it was associated with fewer negative effects on in vitro DM degradability and lower protein binding affinity [164,165]. When the antimethanogenic potential of leucaena extracts with differing average molecular weights was tested, extracts with the highest MW were able engender CH4 to the maximum [161,166]. Table 3 summarizes the data from two different in vitro studies where at the same concentration, the effect of the molecular weight was more pronounced and it had a strong negative correlation with CH4 emissions.

**Table 3.** Influence of molecular weight on methane emissions from two in vitro studies using Kedah-Kelantan cattle rumen fluid.

1 CT: Condensed tannins, 2 MW: Molecular weight of condensed tannins.

Some studies have also analysed the impact of multiple structural features simultaneously from CT extracts. Extracts from multiple sainfoin cultivars and diverse CT sources were analysed for CT structural features such as PD percentage, *cis* flavan-3-ols percentage, and average polymer size (mean degree of polymerization). PD percentage and average polymer size were found play an important role in determining antimethanogenic

potential of CTs, in addition to the actual CT concentration [25,95]. This shows that high reproducibility of the results can be attained by incorporating the structural features of CTs in ruminant nutrition studies.

### *4.3. A Roadmap to Close the Missing Links and Possible Future Directions*

As discussed above, there is an apparent incongruity between measured tannin concentrations and their bioactive effects. This may be explained by a combination of four factors: (a) the variability of tannins and their composition is large both within and amongst species, and it is affected, at least partially, by the environment; (b) most studies have used too few plants and have been conducted under non-controlled environmental conditions or not comparable conditions, to capture the variability of the tannins; (c) many studies have used inadequate or unsuitable analytical techniques (often due to lack of alternatives or resources), which do not capture the structural characteristics; and (d) the studies that investigated the antimethanogenic potential of CTs while accounting for structural attributes are still limited.

To overcome these inconsistencies, in future studies, every aspect that might affect the results, from the growing conditions to growth stage of the plant at harvest, and from sampling to extraction should be carefully considered in future studies. In the absence of laboratory infrastructure for the structural characterisation of CTs, assays to determine their astringency could be employed. Protein precipitation and radial diffusion assays are frequently used to measure the protein binding ability of tannins [167–169]. Furthermore, assays that determine their antioxidative and oxidative (at high pH) behaviour could also be utilized. They have been associated with the negative impact of tannins on the rumen microbial population [90] and their antiherbivore effect [33], respectively. Additional treatments with polyethylene and polyvinylpolypyrrolidone in vitro studies could be employed to elucidate the tannin effect on CH4 emissions, as they bind specifically to tannins. These studies could be instrumental in distinguishing between the effect of tannins and forage chemical composition on CH4 emissions. In vitro fermentation/CH4 production techniques could be useful for screening these forages and determining their adequate supplementation. Using CT extracts of tanniferous forages in in vitro studies can illustrate the structure-activity relationship of CTs with methane emissions more distinctly. Condensed tannin supplementation has been found to impact the diversity and composition of the rumen microbial community [170]. Understanding the dynamics of microbial populations in rumen, and how CT-containing forages influence their abundance and diversity, can provide significant insights into their mode of action. Employing novel techniques such as metagenome and metatranscriptome analysis of the rumen microbiome under CT treatment can help in identifying the microbial population and the functional shifts in rumen microbiome which lead to CH4 abatement [171]. As we gather more information about the relationship between the structure of CTs and their bioactivity, there are prospects for breeding plants with desired concentrations and composition of CTs [105,172]. Molecular approaches have already made it possible for white clover to reach moderate levels of CT in its leaves [173] and efforts are also being made in directions to improve the persistence of TRFs, as may be seen for birdsfoot trefoil [174]. Several questions still remain unanswered, and these are critical for ensuring a comprehensive understanding of the fate of CTs in biological systems.

### **5. Conclusions**

In recent years, there have been remarkable new insights into CT structural diversity and functions with more sensitive analytical methods. However, CT bioactivity is a complex process which results from a multitude of variations occurring simultaneously in plants as well as in their effects in animals. The variability in the results from different studies focuses our attention on the need for developing and adapting a course of action for the investigation of potential of CTs to reduce CH4 emissions. The comparison of CT fingerprints of different species could help us understand not only the factors which define their antimethanogenic potential but also provide a vital framework to assess their interactions with plant constituents and rumen microflora, benefitting overall ruminant health.

**Author Contributions:** S.V. and C.S.M. conceived and conceptualized the manuscript. S.V. took the lead for writing the manuscript, while C.S.M. and F.T. edited for clarity and completeness. All authors provided critical feedback and helped shape the research, analysis and manuscript. Additionally, all authors have read and agreed to the published version of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Deutsche Forschungsgemeinschaft (DFG)—Project number 406534244, grant number MA 8199'/1-1. We acknowledge financial support by DFG within the funding programme Open Access Publizieren.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank Juha-Pekka Salminen for his contributions to improve the quality of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

### *Article* **Technical Efficiency in the European Dairy Industry: Can We Observe Systematic Failures in the Efficiency of Input Use?**

**Lukáš Cechura and Zde ˇnka Ž ˇ áková Kroupová \***

Department of Economics, Faculty of Economics and Management of the Czech University of Life Sciences in Prague, 16500 Prague, Czech Republic<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The paper provides findings on the technical efficiency of the European dairy processing industry, which is one of the most important subsectors of the food processing industry in the European Union (EU). The ability to efficiently use inputs in the production of outputs is a prerequisite for the sustainability and competitiveness of the agri-food sector as well as for food security. Thus, the aim of this paper is to provide a robust estimate of technical efficiency by employing new advances in productivity and efficiency analysis, and to investigate the efficiency of input use in 10 selected European countries. The analysis is based on two-stage stochastic frontier modelling incorporating country-specific input distance function (IDF) estimates and a meta-frontier input distance function estimate, both in specification of the four-component model, which currently represents the most advanced approach to technical efficiency analysis. To provide a robust estimate of these models, the paper employs methods that control for the potential endogeneity of netputs in the multistep estimation procedure. The results, based on the Amadeus dataset, reveal that companies manufacturing dairy products greatly exploited their production possibilities in 2006–2018. The dairy processing industry in the analysed countries cannot generally be characterized by a considerable waste of resources. The potential cost reduction is estimated at 4–8%, evaluated on the country samples mean. The overall technical inefficiency (OTE) is mainly a result of short-term shocks and unsystematic failures. However, the meta-frontier estimates also reveal a certain degree of systematic failure, e.g., permanent managerial failures and structural problems in European dairy processing industry.

**Keywords:** technical efficiency; four-component model; endogeneity; input distance function; metafrontier; stochastic frontier analysis; dairy processing industry; European Union

### **1. Introduction**

A functional, sustainable, competitive, and structurally balanced agri-food sector has an irreplaceable position in the European economy. The importance of this sector in modern society is emphasized by global trends, among which population explosions, migration waves, and climate change can be highlighted. Additionally, the COVID-19 pandemic has strengthened the role of national food security.

The food processing industry is one of the most important sectors of the European economy. According to Eurostat [1], it is the largest manufacturing sector in the European Union (EU), representing 14% of total manufacturing employment, 12% of total manufacturing turnover, and 10% of value added in 2018. The manufacture of dairy products, with a 17% share of total food industry turnover, 12% share of value added, and 10% share of total food industry employment, is one of the most important subsectors of the European food processing industry.

The development of the European food processing industry in recent decades has been affected by several economic, social, and technological trends and challenges, especially the globalization and liberalization of food markets, the global financial crisis, the change

**Citation:** Cechura, L.; Žáková ˇ Kroupová, Z. Technical Efficiency in the European Dairy Industry: Can We Observe Systematic Failures in the Efficiency of Input Use? *Sustainability* **2021**, *13*, 1830. https://doi.org/ 10.3390/su13041830

Academic Editor: Rajeev Bhat Received: 19 January 2021 Accepted: 5 February 2021 Published: 8 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

in consumer preferences towards healthy foods, socially responsible consumption and organic foods, and the implementation of EU regulations focusing on food safety and environmental issues [2,3]. Additionally, the COVID-19 pandemic crisis has created a new era in which the food processing industry is facing various challenges, including the change of consumer purchasing behaviour, transportation network disturbances, labour absenteeism, and the closure of various food manufacturing industries [4]. These trends, connected with a high level of market saturation and concentration of the food retailing sector [5], contribute to a highly competitive environment. To enhance the sustainability and competitiveness of companies/sectors, managers and policymakers have dealt with factors determining productivity, namely, the technology in use, the quality of inputs, the ability to efficiently use inputs in the production of outputs, and the exploitation of economies of scale [6].

These factors have also received special attention in the research of the food processing industry over the last decade. For example, Náglová and Šimpachová Pechrová [7]; Cechura and Hockmann [ ˇ 8]; Rudinskaya [9]; Rezitis and Kalandzi [10]; Špiˇcka [11], Popovi´c and Pani´c [12], and Setiawan et al. [13] investigated the technical efficiency of the food processing industry. Soboh et al. [14] and Dimara et al. [15] dealt with technical and scale efficiency, that was investigated also by Baran [16]. Allendorf and Hirsch [5] analysed technical change and technical efficiency. Cechura et al. [ ˇ 6] assessed the exploitation of economies of scale and production possibilities, along with the impact of technical change. Kapelko et al. [3,17,18], Rudinskaya and Kuzmenko [19], and Cechura [ ˇ 20] investigated productivity dynamics based on technical change, technical inefficiency change, and scale inefficiency change. The majority of these studies are oriented only to one selected country. The exceptions are Kapelko [3], who investigated meat manufacturing firms, fruit and vegetable processing firms, dairy manufacturing firms, and bakery and farinaceous products manufacturing firms in 18 European member states; Allendorf and Hirsch [5], who analysed the dairy and meat processing sectors in eight European countries; Cechura et al. [ ˇ 6], who investigated the milling, fruit and vegetable, dairy, and meat processing sectors in 24 EU member countries; and Soboh et al. [14], who focused on dairy processing firms in six European countries.

With regard to the dairy processing industry, several studies that investigate technical efficiency of European dairy processing industry can be highlighted. Table 1 summarizes these studies. These studies revealed that evaluated on the sample means dairy processing companies in majority of analysed countries greatly exploited their production possibilities in the short-term. However, substantial differences were found between the best and worst producers as well as among countries representing the potential for costs reduction. The studies in Table 1 are in agreement on lower technical efficiency of the Eastern European dairy processing industries compare to the Western European ones. The technical efficiency development was predominantly positive in recent decades and contributed to productivity growth. The need of technical efficiency improvements was prompted by an increase in the price of the main input—milk. Significant short-term changes in technical efficiency were found in the period of the financial crisis. The shortcoming of these studies is that they take into account only transient (short-term) technical efficiency (TTE), which varies across companies because of the shocks associated with new production technologies, human capital, and learning-by-doing. Persistent (long-term) inefficiency (PTE), which could arise due to the presence of rigidity within a firm's organization and production process, is unrecognized in these studies despite the fact that it can lead to the underestimation of technical inefficiency. To the best of our knowledge, there are no empirical studies examining both parts of overall technical efficiency in the European dairy processing industry.

**Table 1.** Recent studies of technical efficiency of European dairy processing industry. Source: authors.

Note: AT denotes Austria, BE denotes Belgium, BG denotes Bulgaria, BIH denotes Bosnia and Herzegovina, CH denotes Switzerland, CZ denotes Czechia, DE denotes Germany, DE denotes Denmark, EE denotes Estonia, ES denotes Spain, FI denotes Finland, FR denotes France, GR denotes Greece, HR denotes Croatia, HU denotes Hungary, IT denotes Italy, LV denotes Latvia, LT denotes Lithuania, LU denotes Luxembourg, PL denotes Poland, PT denotes Portugal, RO denotes Romania, NL denotes the Netherlands, NO denotes Norway, SK denotes Slovakia, SI denotes Slovenia, SR denotes Serbia, SW denotes Sweden, UK denotes the United Kingdom.

The aim of this study is to fill the gap in the empirical literature on the efficiency of the European food processing industry by providing technical efficiency estimates in selected European countries, along with a meta-frontier analysis for cross-country efficiency comparison in the dairy processing industry. In particular, we aim to provide a robust estimate of technical efficiency by employing new advances in productivity and efficiency analysis, and to investigate the efficiency of input use in the analysed countries. The study provides an insight into the decomposition of overall technical efficiency into transient technical efficiency, representing short-term deviations from the frontier, and persistent technical efficiency, which captures systematic deviations from frontier technology.

The study addresses the following research questions: (i) The first relates to the impacts and dynamics of the overall technical efficiency. The aim is to assess whether there is indication that the countries follows a sustainable development path characterized by reduced waste of resource due to inefficient input use. (ii) The second question deals with the sources of overall technical efficiency. The study evaluates whether we can observe systematic failures in the efficiency of input use or whether the deviations from the frontier technologies are due to the transient reasons. Answering these research questions provides important information for policy makers. The structure of the paper is as follows: Section 2 presents the theoretical background of our approach; Section 3 then introduces the model, empirical strategy, and dataset; Section 4 presents the results; and Section 5 discusses our findings and compares them with other studies.

### **2. Theoretical Background**

Technical efficiency was originally defined by Koopmans [23] as the situation where an increase in any output is impossible without a reduction in at least one other output or an increase in at least one input and where a reduction in any input requires an increase in at least one other input or a reduction in at least one output. Debreu [24] and Farrell [25] introduced radial measures of technical inefficiency. Fried et al. [26] suggested an inputconserving orientation. Their measure is defined as the maximum equi-proportionate reduction in all inputs that is feasible with a given technology and outputs. With an output-augmenting orientation, the measure is defined as the maximum equi-proportionate expansion in all outputs that is feasible with a given technology and inputs. A value of unity indicates technical efficiency because no radial adjustment is feasible, and a value different from unity indicates the severity of technical inefficiency.

According to Coelli et al. [27], different methods have been considered for the estimation of the technical inefficiency of the production plan. Two widely used approaches

are data envelopment analysis (DEA), which is non-parametric and deterministic, and stochastic frontier analysis (SFA), which is, on the contrary, parametric and stochastic. There are advantages and disadvantages to each approach (e.g., [28]). The DEA approach is computationally simple, but it is sensitive to outliers since the model ignores measurement error and other sources of statistical noise—all deviations from the frontier are interpreted as the results of technical inefficiency [27].

The stochastic modelling of technical efficiency introduced independently by Meeusen and van den Broeck [29] and Aigner et al. [30] is based on a decomposition of the error term into a symmetric random error and one-sided technical inefficiency. Since the introduction of this frontier modelling within the context of panel-data models, there has been considerable research done in order to extend it to generate consistent and unbiased estimates, see Figure 1. As mentioned by Alem [31], several models have been developed that are based on assumptions about the temporal behaviour of the inefficiency and account for heterogeneity, which can be divided into two categories: observed, which is reflected in the measured variables, and unobserved, which is typically interpreted as the effect of unobservable factors on the outcome of interest [26,32].

**Figure 1.** Development of stochastic frontier models. Source: authors.

Initial studies considering heterogeneity assumed that the time-invariant part of the model represents inefficiency, and the time-varying part represents firm-specific heterogeneity. Later, the firm-specific heterogeneity was assumed to be time-invariant, while the time-variant part was considered as inefficiency [28]. The latest approach emphasises the importance of considering latent heterogeneity in generating an unbiased estimate of time-invariant technical inefficiency, as well as the possibility of efficiency improvement [33,34]. In this approach, the overall technical inefficiency of a producer can be decomposed into transient and persistent parts. Transient technical inefficiency arises as a result of non-systematic managerial failures that can be resolved in the short term. According to Pisulewski and Marzec [35] as well as Njuki and Bravo-Ureta [36], the transient part of inefficiency relates to shocks associated with new production technologies and human capital. Kumbhakar and Lien [37] add that transient technical efficiency can represent the managerial ability to learn by doing. Persistent technical inefficiency represents structural problems in the organization of the production process or a systematic lack of managerial skills [38]. Moreover, it can be an indicator of a non-competitive market condition. Badunenko and Kumbhakar [34] state that persistent inefficiency could not exist in a competitive market, i.e., persistently inefficient firms would not survive in the business. The distinction between persistent and transient technical inefficiency has significant political implications because the persistent part of technical inefficiency is unchangeable without a new policy or fundamental change in the ownership and management of companies. Transient technical inefficiency can be adjusted over time without a major policy change [28].

The properties of the above-mentioned models are described by Kumbhakar et al. [28] and others. A comparison of the majority of SFA panel-data models, proving the sensitivity of technical efficiency estimates to the model specification, is presented by Alem [31] and Rashidghalam et al. [39]. Moreover, Badunenko and Kumbhakar [34] provide a judgement on the reliability of transient and persistent technical efficiency estimates.

### **3. Materials and Methods**

### *3.1. Methodology Used in the Study*

Two-step stochastic frontier modelling is applied to get a robust comparison of the efficiency of input use in the dairy processing industry among the analysed countries. The first step includes the estimation of the country-specific input distance functions, technical efficiency scores, and efficient output levels of dairy processing companies in 10 European member states (Austria (AT), Belgium (BE), Czechia (CZ), Germany (DE), Finland (FI), France (FR), Italy (IT), Spain (ES), Sweden (SW), and the United Kingdom (UK)) using the SFA approach. The second step includes the meta-frontier estimation and a comparison of the technical efficiency among the countries.

### 3.1.1. Input Distance Function

The analysis is based on an assumption that the transformation process is well approximated by an input distance function (IDF) that measures the largest factor of proportionality *ρ* by which the input vector *x* can be scaled down in order to produce a given output vector *y* with the technology existing at a particular time *t* [40]. The IDF is formally defined as

$$D^I(y, \mathbf{x}, t) = \max \left\{ \rho : \frac{\mathbf{x}}{\rho} \in L(y) \right\} \tag{1}$$

where *x* denotes the input vector, *y* denotes the output vector, *t* is a time variable, and *L(y)* represents the input requirement set. For any input-output combination (*x,y*) belonging to the technology set, the input distance function takes a value no smaller than unity. According to Karagiannis et al. [41], a value of unity simply indicates that the input-output combination (*x,y*) belongs to the input isoquant, which represents the minimum input quantities that are necessary to produce a given output vector *y*.

In other words, if *DI* (*y*, *x*, *t*) = 1, the given output vector *y* is produced with the minimum amount of inputs at a given time and with the given technology, and the firm is technically efficient [40]. That is, the IDF provides a measure of technical efficiency since it is reciprocal to the Farrell [25,41] input-based technical efficiency:

$$TE^I = 1/(D^I(y, \mathbf{x}, t))\tag{2}$$

According to Greene [42], the IDF exhibits the following properties: symmetry, monotonicity, linear homogeneity, and concavity in inputs and quasi-concavity in outputs. For the interpretation of the empirical estimates, the duality between the cost and input distance functions is another important property:

$$\mathcal{C}(w, y, t) = \min\_{\mathbf{x}} \left\{ w\mathbf{x} : D^{I}(y, \mathbf{x}, t) \ge 1 \right\} \tag{3}$$

where *w* denotes a vector of input prices. The minimisation problem provides the relation between the derivatives of the IDF and the cost function [43]. In particular, the derivative with respect to the *j*th input gives:

$$\frac{\delta D^I(\mathbf{x}^\*(w, y, t), y, t)}{\delta \mathbf{x}\_j} = \frac{w\_j}{\mathbb{C}(w, y, t)} = r\_j^\*(\mathbf{x}, y, t) \tag{4}$$

That is, the derivative of the input distance function with respect to a particular input is equal to the cost-deflated shadow price of that input. More conveniently, in terms of the log derivative of the distance function, we can rewrite this expression to

$$\frac{\delta \ln D^I(\mathbf{x}^\*(w, y, t), y, t)}{\delta \ln x\_j} = \frac{w\_j \mathbf{x}\_j^\*(w, y, t)}{\mathbf{C}(w, y, t)} = \mathbf{S}\_{j, t} \tag{5}$$

where *Sj,t* is a cost-share of the particular input.

With respect to the output vector *y*, application of the envelope theorem to the minimisation problem *C*(*w*, *y*, *t*) = min *x wx* : *DI* (*y*, *x*, *t*) ≥ 1 leads to

$$\frac{\delta \ln D^I(\mathbf{x}^\*(w, y, t), y, t)}{\delta \ln y\_m} = -\frac{\delta \ln \mathcal{C}(w, y, t)}{\delta \ln y\_m} = e\_{mi, t} \tag{6}$$

Hence, the elasticity of the IDF with respect to the *m*th output is therefore equal to the negative of the cost elasticity of that output, and as such it indicates the importance of output in terms of cost [44].

In this study, we assume that the transformation process can be well approximated by the IDF in a translog functional form. This second-order local approximation of any twicedifferentiable function satisfies Diewert's minimum flexibility requirement for flexible form [35]. The translog input distance function for output (*y*), J-inputs (*x*), and time (*t*) is defined as

 $\ln\text{D}\_{ll}^{I} = \mathfrak{a}\_{0} + \mathfrak{a}\_{m}\ln y\_{il} + 1/2\mathfrak{a}\_{mm}(\ln y\_{il})^2 + \sum\_{j=1}^{I} \gamma\_{mj}\ln y\_{il}\ln\mathfrak{x}\_{j,il} + \sum\_{j=1}^{I} \mathfrak{z}\_{j}\ln\mathfrak{x}\_{j,il} + \mathfrak{z}\_{j}$  $+1/2\sum\_{j=1}^{I}\mathfrak{z}\_{j,k}\ln\mathfrak{x}\_{j,il}\ln\mathfrak{x}\_{k,il} + \delta\_{l}t + 1/2\delta\_{lt}t^{2} + \delta\_{ml}\ln y\_{m,it}t + \sum\_{j=1}^{I} \delta\_{jl}\ln\mathfrak{x}\_{j,it}t$

where subscripts *i,* with *i = 1,2,* . . . *,N*, and *t,* with *t = 1,* . . . *,T*, refer to a certain company and time (year), respectively. *α, β, γ,* and *δ* are vectors of the parameters to be estimated. The symmetry restrictions imply that *βjk = βkj.* The time trend included in the IDF allows for capturing the joint effects of embedded knowledge, technology improvements and learning-by-doing in input quality improvements [45]. Here, *δt* and *δtt* capture the global effect of technical change on the IDF, while *δmt* and *δjt* measure the bias of technical change.

The IDF is homogenous of degree 1 in inputs. According to Sipiläinen [46], it requires

$$\sum\_{j=1}^{J} \beta\_j = 1; \sum\_{j=1}^{J} \beta\_{j\mathbf{k}} = 0; \sum\_{j=1}^{J} \gamma\_{m\mathbf{j}} = 0; \sum\_{j=1}^{J} \delta\_{j\mathbf{t}} = 0. \tag{8}$$

Implying the homogeneity property of the IDF [47], which is imposed by normalising all the inputs by one input, we can rewrite the IDF as

$$\begin{split} & \ln \text{D}\_{\text{It}}{}^{\text{I}} - \ln \text{x}\_{\text{I},\text{it}} = \alpha\_{0} + \alpha\_{m} \ln y\_{\text{it}} + 1/2 \alpha\_{\text{mm}} (\ln y\_{\text{it}})^{2} + \sum\_{j=2}^{\text{I}} \gamma\_{\text{m}j} \ln y\_{\text{it}} \ln \tilde{\text{x}}\_{j,\text{it}} + \sum\_{j=2}^{\text{I}} \beta\_{\text{j}} \ln \tilde{\text{x}}\_{j,\text{it}} + \\ & + 1/2 \sum\_{j=2}^{\text{I}} \sum\_{k=2}^{\text{K}} \beta\_{j\text{k}} \ln \tilde{\text{x}}\_{j,\text{it}} \ln \tilde{\text{x}}\_{k,\text{it}} + \delta\_{\text{l}} \text{t} + 1/2 \delta\_{\text{l}\text{t}} \text{t}^{2} + \delta\_{\text{mt}} \ln y\_{\text{m},\text{it}} + \sum\_{j=2}^{\text{I}} \delta\_{\text{j}\text{t}} \ln \tilde{\text{x}}\_{j,\text{it}} \text{t} \end{split} \tag{9}$$

where ln *x*e*j*,*it* = ln *xj*,*it* − ln *x*1,*it*

Moreover, all variables in logarithm are normalized by their sample mean. In this case, the first-order parameters can be interpreted as output elasticity and input cost shares, evaluated on the sample mean, respectively.

After introducing a statistical error term (*vit*) and latent heterogeneity (*µi* ), and replacing ln *DI it* with inefficiency terms, persistent technical inefficiency (*ηi* ) and transient technical inefficiency (*uit*), that is *ηi* + *uit* = ln *DI it*, the IDF takes the form of a generalized true random effects model (GTRE, [33]):

$$\begin{cases} -\ln \mathbf{x}\_{1,it} = \mathbf{a}\_0 + \mathbf{a}\_m \ln y\_{it} + 1/2\mathbf{a}\_{mm}(\ln y\_{it})^2 + \sum\_{j=2}^l \gamma\_{mj} \ln y\_{it} \ln \tilde{\mathbf{x}}\_{j,it} + \\ + \sum\_{j=2}^l \beta\_j \ln \tilde{\mathbf{x}}\_{j,it} + 1/2\sum\_{j=2}^l \beta\_{j\bar{k}} \ln \tilde{\mathbf{x}}\_{j,it} \ln \tilde{\mathbf{x}}\_{k,it} + \delta\_{\bar{t}}t + 1/2\delta\_{\bar{t}t}t^2 + \\ + \delta\_{\text{mt}} \ln y\_{m,it}t + \sum\_{j=2}^l \delta\_{j\bar{t}} \ln \tilde{\mathbf{x}}\_{j,it}t - \eta\_{\bar{t}} - \mathbf{u}\_{\bar{t}t} + \mu\_{\bar{t}} + \upsilon\_{\bar{t}t} \end{cases} (10)$$

$$\text{where } v\_{\text{il}} \sim N(0, \sigma\_v^2), \mu\_{\text{il}} \sim N^+(0, \sigma\_\text{\textquotedbl}^2), \eta\_{\text{i}} \sim N^+(0, \sigma\_\eta^2), \mu\_{\text{i}} \sim N(0, \sigma\_\mu^2).$$

### 3.1.2. Heterogeneity in Technology

The literature provides broad evidence for significant heterogeneity in dairy processing technology. Since we are estimating a joint country input distance function for the entire food processing industry (due to a data limitation, the low number of observations in an industry does not allow us to estimate country IDF for the dairy processing sector in the majority of countries), we need to consider two types of heterogeneity for processing technology, i.e., the potential existence of inter- and intra-sector heterogeneity. The intrasector heterogeneity is captured by *µi* in (10). To capture the inter-sector heterogeneity, firstorder parameters in (10) are expanded based on dummy variables for four major sectors in the food processing industry (namely the manufacture of dairy products, processing of meat, milling, and manufacture of bakery and farinaceous products):

$$\begin{array}{l} \alpha\_{\rm ms} = \alpha\_{\rm m} + \sum\_{s} \alpha\_{s} d\_{s} \,\forall m\\ \beta\_{\rm fs} = \beta\_{\rm j} + \sum\_{s} \beta\_{s} d\_{s} \,\forall j\\ \delta\_{\rm ls} = \delta\_{\rm l} + \sum\_{s} \delta\_{s} d\_{s} \end{array} \tag{11}$$

where *d* represents dummy variables which account for inter-sectoral differences in technology.

### 3.1.3. Estimation Strategy

Since the endogeneity problem usually frustrates researchers in productivity and efficiency analysis and leads to inconsistent estimates, this study uses methods which control for the potential endogeneity of netputs and thereby allow us to obtain consistent estimates of technology as well as efficiency measures. The study addresses two potential sources of endogeneity (due to the heterogeneity and due to simultaneity of input with technical efficiency) by using the system generalized method of moments (GMM) estimator. In particular, the GMM estimator allows us to deal with the endogeneity bias that arises when one or more explanatory variables are correlated with the error term. This correlation can have different origins, e.g., measurement errors, omitted variables bias, and simultaneity [48].

Moreover, the study compares the GMM estimates with the generalized true random effects model estimates. Analogically to the random effect model, the GTRE assumes that *µi* are independent of explanatory variables. The violation of this assumption can originate from the heterogeneity bias as a kind of omitted variable bias, which is a typical problem of hierarchical data. To deal with this heterogeneity bias, the study applies Mundlak's [49] formulation and adds group-means for each time-varying explanatory variable in the first-order level.

Since the GMM model deals with both sources of endogeneity, it represents our model of first choice. If the GMM model does not provide consistent estimates, the GTRE model with Mundlak's adjustment is our second choice. The standard GTRE model allows for an overall model comparison.

The estimation of the GTRE model and the GTRE model with Mundlak's adjustment is undertaken as a multistep procedure. We follow Kumbhakar et al. [50] and rewrite the model in (10) as

$$\begin{cases} -\ln \mathbf{x}\_{1,it} = \mathbf{a}\_0^\* + \mathbf{a}\_m \ln y\_{it} + 1/2\mathbf{a}\_{mm}(\ln y\_{it})^2 + \sum\_{j=2}^l \gamma\_{mj} \ln y\_{it} \ln \tilde{\mathbf{x}}\_{j,it} + \\ + \sum\_{j=2}^l \beta\_j \ln \tilde{\mathbf{x}}\_{j,it} + 1/2\sum\_{j=2}^l \beta\_{jk} \ln \tilde{\mathbf{x}}\_{j,it} \ln \tilde{\mathbf{x}}\_{j,it} + \delta\_{lt}t + 1/2\delta\_{lt}t^2 + \\ + \delta\_{mt} \ln y\_{m,lt}t + \sum\_{j=2}^l \delta\_{jt} \ln \tilde{\mathbf{x}}\_{j,it}t + \alpha\_i + \varepsilon\_{it} \end{cases} \tag{12}$$

where *α* ∗ 0 = *α*0 − *E*(*ηi*) − *E*(*uit*), *αi* = *µi* − (*ηi*) − *E*(*ηi*)) and *εit* = *vit* − (*uii* − *E*(*uit*)).

This specification ensures that *αi* and *εit* have zero mean and constant variance. The multistep procedure consists of three steps. In step 1, standard random effect panel regression is used to estimate *β*, *γ, δ*, *αm*, and theoretical values of *αi* and *εit*, denoted by *α*ˆ*i* and *ε*ˆ*it*. In step 2, the transient technical inefficiency, *uit*, is estimated using *ε*ˆ*it* and the standard stochastic frontier technique with assumptions: *vit* ∼ *N*(0, *σ* 2 *v* ), *uit* ∼ *N*+(0, *σ* 2 *u* ). In step 3, the persistent technical inefficiency, *ηi ,* is estimated using and the stochastic frontier model with the following assumptions: *ηi* ∼ *N*+(0, *σ* 2 *η* ), *µi* ∼ *N*(0, *σ* 2 *µ* ). These steps are done in the SW STATA 14.0.

Furthermore, the total technical efficiency (OTE) is quantified based on Kumbhakar et al. [50]:

$$OTE\_{it} = \exp(-\hat{\eta}\_i) \times \exp(-\hat{u}\_{it}).\tag{13}$$

The GMM model extends this estimation procedure. The study follows the four-step procedure of Bokusheva and Cechura [ ˇ 51]. In step 1, the two-step system generalized method of moments estimator [52,53] is used to obtain consistent estimates of the IDF parameters. The Arellano and Bover [52]/Blundell and Bond [53] approach builds a system of two equations, the original equation (in levels) and the transformed one (in differences), and employs two types of instruments: level instruments for the differenced equations and lagged differences for the equations in levels [54]. The validity of these instruments is tested by the Hansen J-test [55], which analyses the joint validity of the instruments, and the Arellano-Bond test [56], which analyses the autocorrelation in the idiosyncratic disturbance term (*vit*) that could render some lags invalid as instruments. In step 2, residuals are used from the system GMM level equations to estimate a random effects panel model employing the generalized least squares (GLS) estimator. The transient and persistent technical inefficiency is estimated in steps 3 and 4 based on the same procedure as described above. These estimates are also done in the SW STATA 14.0.

### *3.2. Data Used in the Study*

The study uses an unbalanced panel data set drawn from the Amadeus database collected by Bureau van Dijk—Moody's Analytics company. The database contains information on around 21 million European companies and provides detailed information about company financials in a standard format, financial strength indicators, sectoral activities, and corporate structures. Moreover, the database is unified between different countries to guarantee the comparability of data. The scope of the database and the comparability of data are the reasons why this database is widely used in economic research. These are examples of studies of technical efficiency of dairy processing industry based on Amadaus database: Kapelko and Oude Lansink [21], Cechura et al. [ ˇ 6], Soboh et al. [14], and Špiˇcka [11].

This study uses information from the final accounts of companies whose main activity was food processing (Division C10: Manufacture of food products according to the Statistical classification of economic activities in the European Community, abbreviated as NACE) in the period from 2006 till 2018. Moreover, the study is concentrated on an analysis of C10.5: Manufacture of dairy products (including subgroups: 10.51, 10.52). The study focuses on 10 countries which together represents 76% of European food processing turnover according to Eurostat [1]. Specifically, the analysis uses the following data: output (y), labour (xL), capital (xC), and material (xM). The output is represented by operating revenue (turnover) deflated by the sectoral index of food processing prices (EU level or country level if it was disposable; 2010 = 100) and changes in a company's stock. Labour is represented by the cost of employees deflated by the index of producer prices in the industry (country level; 2010 = 100). Capital is the book value of fixed assets deflated by the index of producer prices in the industry (country level; 2010 = 100). Finally, material is the total cost of materials and energy deflated by the index of producer prices in the industry (country level; 2010 = 100). The source of price indexes is the EUROSTAT database.

Since not all information can be found in the database, only those companies with non-zero and positive values are used for the variable of interest. Moreover, companies with less than three consecutive observations are rejected from the dataset. This procedure decreases the problem associated with the entry and exit of producers from the database

and allows the use of the GMM estimator with a sufficient number of lagged instruments. The study is constrained by that to the use of an unbalanced panel data set containing 11,605 food processing companies with 95,003 observations in the first step: country-specific IDF estimation. The second step: meta-frontier IDF estimation uses only data of dairy processing companies, and the structure of this data sub-set is presented in Table 2.

Note: I is the number of dairy processing firms (10.5), N is the total number of observations (10.5) RS1 is the revenue share (in %) of Amadeus sample food processing firms on the total revenue of food processing firms, RS2 is the revenue share (in %) of Amadeus dairy processing firms on the total revenue of dairy processing industry; AT denotes Austria, BE denotes Belgium, CZ denotes Czechia, DE denotes Germany, ES denotes Spain, FI denotes Finland, FR denotes France, IT denotes Italy, SW denotes Sweden, UK denotes the United Kingdom.

### **4. Results**

The country-specific IDF estimates (country-specific IDF estimates are available on request from the authors or in [57]) reveal that the technology of the food processing industry exhibits common patterns in the analysed European countries. Common features are found in the case of cost elasticity, as well as in the cost structure. As expected, and in line with the information provided by the dataset, the highest cost shares have material inputs between 0.50 and 0.64. The labour cost share is approximately 30%, with the lowest cost share in Germany (25.8%) and the highest in France (38.3%). The rest is represented by the share of capital. Moreover, food processors in general cannot benefit greatly from the exploitation of economies of scale since the size seems to be optimal in the majority of cases. The elasticities of output, which reveal whether the technology exhibits increasing, constant, or decreasing returns to scale (i.e., whether a proportionate increase in all inputs results in a larger, equal, or less than proportionate increase in the aggregate output, respectively [58]), are close to one in the majority of countries, evaluated on the sample mean. The deflection towards a higher absolute value than one is a characteristic of the Nordic countries (Finland and Sweden).

As far as heterogeneity is concerned, the results do not reveal any significant heterogeneity effect on output elasticities for the majority of countries. In other words, inter-sectoral heterogeneity in cost elasticity is not a common feature of European food processing sectors. When the focus is primarily on the manufacture of dairy products, an exception can be found in Spain, where an absolute value of output elasticity slightly higher than one reveals that dairy processors exhibit moderate economies of scale, evaluated on the country sample mean. Inter-sectoral heterogeneity is not a statistically significant feature of food processing industries, even in the case of cost shares. This fact holds true especially for labour. The heterogeneity parameters of labour are statistically significant (at the 10% level) only in the Spanish and Swedish dairy industries (without a common pattern). More frequent occurrence of significant inter-sectoral heterogeneity is revealed in the case of material cost shares. The dairy processing sector is characterized by higher material cost shares in the majority of countries compared to the rest of the food processing sectors.

The common patterns of the food processing industry in the analysed European countries are also revealed in the term of technical change. The technology exhibits positive technical change in most of the analysed countries, and this is accelerating over time. With a focus on the dairy processing industry and an evaluation on the country sample means, a technological regression is revealed only in Belgium and Finland. This may indicate a low level of innovation (especially process one) in these countries and may imply significant structural weakness [59]. In general, the lack of innovation activities may result in lower competitiveness. On the contrary, a successful innovation process resulting in cost diminution is an important source of competitive advantage. These results have important

political implications. Strengthening the innovative activity of companies is unlikely to be possible without investment incentives and continued support for the building of an innovation-friendly business environment. Supporting the innovation activity of food processors is an important source not only of profitability and competitiveness but also of food security, safety, and sustainability.

Table 3 summarises parameter estimates of the stochastic meta-frontier input distance function models in the three alternative specifications for the dairy processing industry in 10 European countries. Considering theoretical consistency, the IDF estimates should satisfy the properties of monotonicity and concavity in inputs. The violation of these assumptions may result in misleading values of computed elasticities and technical efficiency scores. The monotonicity assumption is met if the IDF is non-increasing in output and non-decreasing in inputs. In order to verify this property, it is sufficient that the first-order parameters are negative for the output (*βm* ≤ 0) and positive for the inputs (*βj* ≥ 0 for *j = L,K*) and *βL + βM* < 1, where *L* stands for labour and *M* represents material [60]. Concavity in inputs requires that the Hessian matrix of second-order derivatives of the IDF of the function with respect to the inputs is negatively semidefinite, according to Diewert and Wales [61]. This is fulfilled on the sample mean, if *βjj + β 2 j –βj* ≤ 0 *j = L,M.* Table 3 proves that these conditions hold for all model specifications, evaluated on the sample mean.

Note: GMM model–AR (2) = −0.06 and Hanset test = 347.25 (No. of instruments–226).

The first-order parameters of the IDFs are highly significant even at the 1% significance level and have a similar pattern, indicating the high material intensity of the dairy industry. The share of material in the total input is about 60%, the share of labour is about 30%, and the output elasticity varies between (−0.98) and (−0.89). The time parameter (*δt*) has a low negative value representing the cost decrease with a decelerating rate (*δtt >* 0) over the analysed period. However, the second-order time parameter (*δtt*) is statistically significant only at the 10% significance level in the GMM model. The parameters of biased technical change are not statistically significant even at the 10% significance level. The rest of the second-order parameters have a similar pattern as well.

The results reveal high overall technical efficiency. Evaluated at the sample averages, the overall technical efficiency is between 92% and 96%, that is, dairy processors can reduce

their costs from 4% up to 8%, subject to the model specification, to produce the same volumes of outputs. As far as the standard deviations are concerned, we may observe that the dairy producers operate very close to the estimated mean value of technical efficiency. Figure 2 shows that the models provide a similar overall technical efficiency distribution.

**Figure 2.** Kernel density comparison (overall technical efficiency). Source: authors' own calculation.

Persistent technical inefficiency is estimated only in the GMM model and is quite low. That is, with similar estimates of transient technical efficiency, we do not see significant differences in overall technical efficiency among the model specifications. These results suggest that the efficiency estimates seem to be robust.

Table 4 provides the statistical characteristics of the country-specific overall technical efficiency and its decomposition into transient and persistent parts. The overall technical efficiency estimates indicate that, evaluated at the country sample means, dairy processing companies can reduce their costs from 5% up to 24% while producing the same volumes of output.

As far as intra-sectoral differences are concerned, the overall technical efficiency is relatively dense around the mean and the distribution is skewed to higher values. In the majority of the analysed countries, the interquartile ranges are quite narrow, see Figure 3. The outliers in boxplots indicate high polarization between the most technically efficient producer and the least technically efficient one. The highest polarization is evident in Austria. On the other hand, the lowest polarization can be seen in Czechia and Finland. As the overall technical efficiency is associated with cost savings, it can be assumed that Austria, with the highest polarization, will face structural changes in the future. That is, we can expect changes in business activities or even the cessation of production for firms with the lowest technical efficiency and resource reallocation to more successful producers.

**Figure 3.** Country-specific overall technical efficiency (GMM). Source: authors' own calculation. Note: AT denotes Austria, BE denotes Belgium, CZ denotes Czechia, DE denotes Germany, ES denotes Spain, FI denotes Finland, FR denotes France, IT denotes Italy, SW denotes Sweden, UK denotes the United Kingdom.

The overall technical efficiency consists of the transient and persistent technical efficiency. The persistent part of overall technical efficiency indicates systematic failures or structural problems and unsuitable factor allocations, which are difficult to change over time, as well as non-competitive market conditions. However, the persistent technical efficiency is statistically significant only in the cases of Austria and Italy. Table 4 shows that the average persistent technical efficiency scores are about 91% in these two countries. A slightly higher variability for persistent technical efficiency is revealed in Austria. However, the Italian dairy sector is characterized by the highest polarization between the most efficient producer and the least efficient. This suggests that many firms systematically fail to catch up to best practices and thus show considerable resource inefficiency as compared to firms operating on the technological frontier.

Table 4 also provides country-specific transient technical efficiency scores. The intrasectoral differences of transient technical efficiency are pronounced especially in Austria, Italy, and Germany. However, the highest polarization is revealed in France, followed by Germany, where the developments of transient technical efficiency over time have a similar pattern and the lowest values are connected with the 2009 crisis; see Figure 4. The fluctuations of technical efficiency in the short-term may also be the result of shocks associated with the introduction of new technologies or changes in human capital. Considering the development presented in Figure 4, we can conclude that the most important feature is an increasing trend in transient technical efficiency in the analysed period.

**Figure 4.** Country-specific transient technical efficiency development (GMM). Source: authors' own calculation. Note: AT denotes Austria, BE denotes Belgium, CZ denotes Czechia, DE denotes Germany, ES denotes Spain, FI denotes Finland, FR denotes France, IT denotes Italy, SW denotes Sweden, UK denotes the United Kingdom.

Table 5 presents meta-frontier technical efficiency estimates and indicates that the overall technical efficiency is high in all countries and there are no considerable differences among them. The lowest value of overall technical efficiency can be observed in Sweden, evaluated on mean values. Moreover, the Swedish dairy industry can be characterized by the highest variability measured by standard deviation. On the contrary, the highest average value of overall technical efficiency is found in France and Finland.

**Table 5.** Meta-frontier overall technical efficiency and its decomposition (GMM). Source: authors' own calculation.

Figure 5 illustrates the competitive position of selected countries in terms of a countryaverage technical efficiency comparison with the average value of the whole set. As can be seen, the Swedish, Belgian, and German dairy processing industries are failing to catch up with the best-practice technology. While in Belgium and Germany the dairy processing industry lags behind in both transient (TTE) and persistent (PTE) technical efficiency, in Sweden the transient technical inefficiency poses a greater problem for the dairy processing industry than the persistent component.

**Figure 5.** Meta-frontier technical efficiency comparison (GMM). Source: authors' own calculation. Note: AT denotes Austria, BE denotes Belgium, CZ denotes Czechia, DE denotes Germany, ES denotes Spain, FI denotes Finland, FR denotes France, IT denotes Italy, SW denotes Sweden, UK denotes the United Kingdom.

The loss of resources due to structural problems and permanent managerial failures in the production process is pronounced in all countries and reaches a similar level. However, it holds in all cases that the persistent inefficiencies are not large, evaluated on the sample means. On the other hand, the outliers in boxplots (Figure 6) indicate high polarization between the most long-term technically efficient producer and the least technically efficient one in Belgium, Germany, and the United Kingdom. In other words, it reveals the existence of firms that are systematically failing to catch up to best practices and thus show considerable resource inefficiency. Different paths can be followed to eliminate this waste of resources and promote the sustainability of the dairy processing industry. According to Pisulewski and Marzec [35], management should focus on changes in the organization of the production process and managerial competency. According to Dimara et al. [15], policy makers should provide access to biotech innovations and food supply networks.

**Figure 6.** Meta-frontier persistent technical efficiency comparison (GMM). Source: authors' own calculation. Note: Reference line represents the sample average transient technical efficiency score; AT denotes Austria, BE denotes Belgium, CZ denotes Czechia, DE denotes Germany, ES denotes Spain, FI denotes Finland, FR denotes France, IT denotes Italy, SW denotes Sweden, UK denotes the United Kingdom.

### **5. Discussion**

The country as well as meta-frontier technical efficiency estimates revealed high overall technical efficiency for all countries. The only exception is the country estimate for Austria. These results are consistent with the literature devoted to technical efficiency analysis in the food processing industry. For example, Cechura et al. [ ˇ 6] found high technical efficiency in the food processing industry in the analysis of entire EU food processing sectors. In addition, their meta-frontier analysis provides similar results as well. The authors show that the differences in average technical efficiency are not large for all EU member countries, including Serbia, even when there are huge differences between the best and worst food processors. High technical efficiency of the top 10% of food processors is a common feature for all countries in all analysed sectors. Similar findings were provided for Czechia by Cechura and Hockmann [ ˇ 8] and Špiˇcka [11]. On the other hand, lower technical efficiency scores were estimated for the Czech food processing industry by Náglová and Šimpachová Petrová [7]. The differences are due to the use of both a different data set and a different model specification. The authors analysed the manufacturing of food products and beverages. Moreover, they did not consider heterogeneity among sectors, which may result in higher inefficiency scores. The treatment of the unobserved heterogeneity component in the model specification is a potential source of the slightly lower technical efficiency scores obtained by Soboh et al. [14] as well.

Our technical efficiency complements the research of technical efficiency in the food processing industry by the decomposition of overall technical efficiency into transient and persistent components. Transient technical efficiency dominates in the technical efficiency country estimates. That is, the overall technical efficiency estimated for each country is largely due to the transient technical efficiency. Persistent technical efficiency is pronounced

only in Austria and Italy. These findings suggest that in the majority of countries, the inefficiencies in input use do not have a systematic character. This holds true for the food processing sector as a whole, since the country estimates are based on a joint IDF function for the four processing sectors. That is, the estimated inefficiencies have non-systematic sources such as non-systematic managerial failures, market shocks, shocks associated with new production technologies, and human capital. However, the meta-frontier estimates that compare the dairy frontier technologies among the countries show a certain degree of systematic failure. In particular, the average persistent technical efficiency is high and similar for all countries. The same holds true for transient technical efficiency. Thus, the meta-frontier analysis shows high overall technical efficiency in the dairy processing sectors in all analysed countries. These findings are fully in line with the meta-frontier analysis provided by Cechura et al. [ ˇ 6].

### **6. Conclusions**

The study aimed to (i) assess whether there is indication that the countries follows a sustainable development path characterized by reduced waste of resources due to inefficient input use and (ii) whether we can observe systematic failures in the efficiency of input use or whether the deviations from the frontier technologies are due to the transient reasons. The results revealed that both components of overall technical efficiency are high. That is, we cannot observe considerable systematic failures in efficiency of input used as well as inefficiencies due to the transient reason, evaluated on sample means. This suggests that the European dairy processing industry as a whole seems to be highly competitive and the companies highly efficient. On the other hand, the figures show that there are companies that are falling behind and may leave the market.

The findings of this study seem to be of high relevance for policy makers since there is still limited knowledge about the inefficiency of input use in the processing industry and their sources and thus the potential contribution the sustainability improvements in the whole food value chain.

There are some limitations of the conducted research. This is especially related to the fact that the Amadeus data allows to work only with aggregate output and thus the diversification or economies of scope cannot be considered. On the other hand, the advantages are the employment of new advances in productivity and efficiency analysis and thus the robust estimate of technical efficiency in European dairy industry.

Future research should focus on investigation of the socioeconomic and environmental factors that can explain differences in persistent as well as transient inefficiency between dairy processing companies. Understanding these factors has important political implications especially in food processing industry that face the challenge of ensuring a sufficient food supply from limited resources for a growing population. Industrial policy agents can use the knowledge of these factors to make a business environment that improves sustainability and competitiveness of food processing industry.

**Author Contributions:** Conceptualization, L.C. and Z.Ž.K.; methodology, L. ˇ C. and Z.Ž.K.; validation, ˇ L.C.; formal analysis, L. ˇ C. and Z.Ž.K.; investigation, L. ˇ C. and Z.Ž.K.; resources, L. ˇ C. and Z.Ž.K.; ˇ writing, L.C. and Z.Ž.K.; visualization, Z.Ž.K.; supervision, L. ˇ C.; project administration, L. ˇ C.; funding ˇ acquisition, L.C. Both authors have read and agreed to the published version of the manuscript. ˇ

**Funding:** The VALUMICS project "Understanding Food Value Chain and Network Dynamics" received funding from the European Union's Horizon 2020 research and innovation programme, under grant agreement No. 727243. https://valumics.eu/.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data was obtained from the CULS and are not available from the authors. The data can be bought from Bureau van Dijk—Moody's Analytics company.

**Acknowledgments:** The results are the part of the solution of the VALUMICS project "Understanding Food Value Chain and Network Dynamics". The project received funding from the European Union's Horizon 2020 research and innovation programme, under grant agreement No. 727243. https://valumics.eu/.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

### **References**

### *Article* **Environmental Impacts of Milking Cows in Latvia**

**Janis Brizga 1,\* , Sirpa Kurppa 2 and Hannele Heusala 2**

**Abstract:** Increasing pressures surrounding efficiency and sustainability are key global drivers in dairy farm management strategies. However, for numerous resource-based, social, and economic reasons sustainable intensification strategies are herd-size dependent. In this study, we investigated the environmental impacts of Latvia's dairy farms with different management practices. The herd size-dependent management groups varied from extensively managed small herds with 1–9 cows, extending to stepwise more intensively managed herds with 10–50, 51–100, 100–200, and over 200 milking cows. The aim is to compare the environmental impacts of different size-based production strategies on Latvia's dairy farms. The results show that the gross greenhouse gas emissions differ by 29%: from 1.09 kg CO2 equivalents (CO2e) per kg of raw milk for the farms with 51–100 cows, down to 0.84 kg CO2e/kg milk for farms with more than 200 cows. However, the land use differs even more—the largest farms use 2.25 times less land per kg of milk than the smallest farms. Global warming potential, marine eutrophication, terrestrial acidification, and ecotoxicity were highest for the mid-sized farms. If current domestic, farm-based protein feeds were to be substituted with imported soy feed (one of the most popular high-protein feeds) the environmental impacts of Latvian dairy production would significantly increase, e.g., land use would increase by 18% and the global warming potential by 43%. Environmental policy approaches for steering the farms should consider the overall effects of operation size on environmental quality, in order to support the best practices for each farm type and steer systematic change in the country. The limitations of this study are linked to national data availability (e.g., national data on feed production, heifer breeding, differences among farms regards soil type, manure management, the proximity to marine or aquatic habitats) and methodological shortcomings (e.g., excluding emissions of carbon sequestration, the use of proxy allocation, and excluding social and biodiversity impacts in life-cycle assessment). Further research is needed to improve the data quality, the allocation method, and provide farm-size-specific information on outputs, heifer breeding, manure storage, and handling.

**Keywords:** LCA; dairy; environmental impacts; intensification

### **1. Introduction**

Dairy production systems are important and complex sources of environmental impacts including global warming, nutrient losses, water and land use [1]. The sustainable intensification of animal production is generally promoted as a way to mitigate these environmental impacts while increasing productivity [2,3]. Nevertheless, the dairy sector has shown a global trend toward intensification [4].

During the last decades following the global trend, dairy cow management systems in Latvia are intensifying rapidly [5–7], i.e., the role of grazing has decreased, and indoor feeding and high protein diets have increased. This intensification process has been in parallel with the concentration of the dairy industry, where small farms leave the market whereas the most competitive farms continue growing and increasing their herds [7]. Nevertheless, the dairy systems in Latvia are still very diverse and are largely determined by the number of cows in a herd. Eighty-four percent of the dairy farms are small—less

**Citation:** Brizga, J.; Kurppa, S.; Heusala, H. Environmental Impacts of Milking Cows in Latvia. *Sustainability* **2021**, *13*, 784. https://doi.org/10.3390/su13020784

Received: 26 October 2020 Accepted: 11 January 2021 Published: 15 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

than 10 animals per farm—producing 18% of the country's raw milk. Most of them are self-sufficient farms that do not deliver dairy products to the market. A similar situation exists in many Eastern European countries as well as in developing countries [8–10].

Data obtained from a survey in Latvia shows significant differences in milk yield depending on the size of the herd (see Table 1). Milk yield depends on diet and the cows' ability to consume dry matter (DM) [11]. While consumption of DM is affected by the quality and digestibility of feed and the optimum ratio of roughage and feed concentrate in the feed ration [12], diets, on the other hand, are linked to the size of the herd affecting the milk yield. The average milk yield of the farms that graze their dairy cows is lower than that of the farms who do not graze their dairy cows [13].

**Table 1.** Description of the technical management routes.

Apart from the milk yields and diets of different size dairy farms, they also differ by other parameters, e.g., the milking technologies and cowsheds (heat-insulated and uninsulated) used [13], as well as the cow breeds. Two dairy cow breeds account for almost 75% of the total dairy cows in Latvia—Latvian Brown, and Black and White Holstein. However, there is no reliable data that shows how these factors create differences in a farm's output among different farm sizes.

Empirical research on the environmental effects of intensification shows that increased dairy production per ha leads to increased impacts per ha, but environmental impacts per kg of raw milk are not as clear [15]. Other studies show that increasing milk production per ha may lead to reduced global warming potential (GWP) per kg of milk [16,17] and improvements in feed efficiency, fertility, and cow longevity are important parameters to reach increased milk yield and lower climate and land footprints per kg of milk [18,19].

This study aims to compare the environmental impacts of raw milk production by different size dairy farms in Latvia. The goal was to use the life-cycle impact assessment (LCA) method to gain insights into differences among farms in terms of their environmental impacts. In addition, the effect of protein feed substitution is also assessed. These results will help better tailor environmental policies to different size farms, and thus decrease the environmental impacts of dairy farming.

### **2. Materials and Methods**

Life-cycle assessment (LCA) is a widely acknowledged and standardized method to evaluate environmental impacts during the entire life cycle of a product (ISO, 2006). Therefore, we developed the LCA model to capture the most important interactions in complex dairy production systems and determine cradle-to-farm gate (all processes involved in milking farms and production of inputs used in farms) environmental impacts of milking cows in Latvia.

Thus, the model describes the diversity in diets and the main inputs (energy, water, and minerals) for different milking farm management strategies in Latvia. More specifically, the system boundaries of this study include the following processes and flows: emissions from manure storage and application and enteric fermentation by ruminants, direct and indirect emissions from animal feed, energy, and water use. Diet composition is an input to this module which includes the cultivation and processing emissions of all feed resources, such as grass (including grazing and hay), alfalfa and maize silage, fodder, roots, as well as molasses, rapeseed expeller and meal, salt, minerals, and vitamins. However, infrastructure and machinery are excluded from the scope of this study.

Although heifers are an important part of the dairy production system, the growing of heifers could not be included in this study as the database used does not have information on heifer rearing processes. Heifer rearing would need a specific study because there are also big differences in the existing production practices. Small farms mostly grow the heifers in the same farm, heifer farming being a side role of the farm. Big farms, however, mostly buy heifers externally from farms that specifically concentrate on heifer farming.

In LCA, the environmental impact is expressed per functional unit, which in this study was one kg of fat and protein corrected milk (FPCM) leaving the farm-gate, which allows us to compare the environmental impacts of different farming systems—the farms with different herd sizes. The following equation was used to calculate FPCM [20]:

### FPCM (kg/a) = milk (kg/a) × ((0.1226 × fat %) + (0.0776 × protein %) + 0.2534)) (1)

Dairy farms generate several outputs to the technosphere. The herd at the average Latvian dairy farm consists of nine dairy cows. Dairy cows which are severely injured or old are slaughtered and hardly any male animals are kept, while most female calves are kept and raised for herd replacement. Most of the male calves and a small number of the female calves which are not needed for herd replacement are sold shortly after birth. Therefore, the environmental impacts in this study have been allocated across these co-products utilizing a proxy mass-based allocation approach—raw milk (85.95%), cows for slaughter (12.35%), and calves (1.7%) [21]. In this study these percentual allocations are applied equally across all the farms as we are missing farm-size-specific information.

Most of the data were collected from secondary sources (see Table 1). Data on herd sizes were obtained from the Central Statistics Bureau (CSB) database [6]. Information about the protein and fat content in the milk for different size farms, used to calculate the FPCM, was obtained from the national Agriculture Data Center (ADC) [14]. These are reliable data sources that are also used in the national inventory to calculate agricultural gross greenhouse gas (GHG) emissions [22].

Diet composition as kg dry matter per cow (see Table 2) and average daily nutrient provision in feed rations for cows in different size farms (see Table 3) were obtained from the survey done by Degola et al. [23], where they obtained data about feed intake values in different size herds from 24 animal feeding experts from the Latvian Rural Advisory and Training Centre. It should be mentioned that the dietary intakes in Table 2 represent average data for each of the fodder in a particular farm group. In smaller farms, there is a wide diversity in feeding strategies, and thus the daily intake in DM differs from the daily DM consumption per cow in Table 3. However, feed related impact calculations are based on the total daily DM consumption from Table 3. Emission factors of feed production were obtained from a life-cycle inventory data from the EcoInvent 3 database [21].

**Table 2.** Average daily feed intake per cow in different management groups (kg DM).

Source: Degola, Cielava, Trupa and Aplocin ¯ , a [23].

**Table 3.** The average daily nutrient provision in feed rations for cows in different management groups.

Source: Degola, Cielava, Trupa and Aplocin ¯ , a [23].

We also had to use several assumptions to differentiate the farms. For heating energy use we assumed that the smallest farms (1–9 cows) were using non-heated barns, while other farms (>9 cows) were using wood biomass for the milking room heating [24]. Electricity consumption was differentiated based on the assumption that smaller farms are not using power for ventilation, milk cooling, and milking (see Table 4). However, we have to acknowledge that some of the information is missing and not covered by this study, e.g., the proximity of the farms to marine or aquatic habitats or a difference in soil quality, fertility, type, and topography among farm sizes.

The cumulative life-cycle inventory data were assessed at the mid-point level using the ReCiPe 2016 (H) V1.02 method. Following the most relevant environmental impact categories for dairy production proposed by the European Dairy Association [25], assessments were made of:

Additionally, we analyzed freshwater ecotoxicity as we believe it is one of the key issues for sustainable agriculture and that dairy farmers can perform actions to reduce the toxic impact of dairy production on ecosystems [11]. Sometimes insecticides are also used directly on the animal to protect them from insects during hot summer days. However, because of the lack of reliable data, this is not included in the study.

With regards to the GWP, in this study, we tried to cover all the major emissions, but because of the lack of reliable data, soil carbon sequestration is not included in the calculation of the carbon footprint of crop production in this study. Similarly, nitrogen fixation by plants is excluded from the study as it is a soil quality dependent process and soil data for the different farm types is not available.

Annual methane emissions from enteric fermentation per cow were calculated using the following equation, provided by the United Nations Food and Agriculture Organization (FAO) [2]:

$$\text{Annual CH}\_4 \text{ emission} = \text{((DMa \times Ym/100))} \times \text{(18.55/55.65)}\tag{2}$$

where:

Methane emissions from enteric fermentation were calculated to be from 85.8 in the small farms to 148 kg CH4 cow−1 year−1 in the largest farms (>200 cows). The average annual national methane emission factor was calculated to be at 118 kg CH4 cow−1 , while the national methane emission factor used for the year 2016 in the National Inventory was 136 kg CH4 cow−1 year−1 [22]. Latvia uses a higher emission factor for dairy cows based on a different feeding situation that is not fully characterized as stall-fed (set for Tier 1). Besides, the digestibility used for calculations of emission coefficient is lower (65% against 70% for Tier 1).

Methane emissions depend on the composition and amount of manure produced (anaerobic conditions, found at the bottom of deep lagoons, increase methane emissions) and the type of storage used [2]. Methane productivity is also determined by the quantity of volatile solids in the manure which depends on the digestibility of feed—low digestibility feed results in a high amount of volatile solids in manure which are then converted to methane. In Latvian small farms, milk cows are still mainly stanchioned, producing farmyard manure, but with increasing intensification there is a gradual transition to the producing of liquid manure [27]. It is expected that the diets of cows (high levels of grass and silage) on the smaller farms will result in higher N excretion because the amount and form of the N excreted is affected by the concentration and ruminal degradability of crude protein in the diet. However, it should be noted that these factors are not fully considered

in this study as there are no data available on manure management for the different sizes of farms in Latvia. Therefore in this study, annual methane emissions of manure management and N excretion rates for dairy cows were taken from the national inventory and are set as 17.28 kg CH4 cow−1 year−1 and 113.9 kg N cow−1 yr−1 respectively [22].

### **3. Results**

The analysis of the environmental impacts of milking cows by farms practicing different management systems in Latvia shows that the relationship between farm management practices/size and environmental impacts is not straightforward, i.e., in some impact categories (water footprint and land use) the small farms create the most burden, but in other categories (e.g., global warming, marine eutrophication, fossil resource scarcity, terrestrial acidification) the highest impacts are caused by the mid-sized farms (see Table 5).

**Table 5.** Results of the characterization: comparison between environmental pressures of different size dairy farms in Latvia per kg of FPCM.

For all the farms, the largest GWP is associated with enteric fermentation and manure (see Figure 1). For smaller farms, these factors account for 61% of the GWP, but for farms with 200 or more cows these factors account for 68% of all the impact per kg of raw milk produced. The second most important process is alfalfa-grass silage, which peaks (15% of GWP) in the mid-sized farms. The next most important factor is concentrates (in our case pea and oat), which combined accounted for 9 to 15% of GWP. These feed-related GHG emissions are mostly from the energy used in field management and harvest. Emissions associated with electricity production are responsible for 2 to 3% of GWP, but in this case, smaller farms have relatively smaller electricity consumption than larger farms. Results for GWP depend mainly on different diet efficiencies and energy use. Global warming potential in our study varies from 1.1 kg CO2e/kg of milk from the farms with 51–100 cows to 0.84 kg CO2e/kg of milk from farms above 200 cows. Thus, the national average gross GWP is 0.9 kg CO2e/kg of raw milk, but GWP per cow is between 5200 kg of CO2e in the smallest farms and 7700 kg of CO2e in the farms with 51 to 100 cows (see Table 6). The total cradle-to-farm gate GWP of the milking cows in Latvia is estimated to be 1 million tons of CO2e.

**Figure 1.** Global warming potential (t CO2e/t milk) for different dairy systems in Latvia.

**Table 6.** Annual GWP and land use for different management groups.

One of the important impact categories is land use, which from all the parameters analyzed has the highest coefficient of variation among the different farm sizes (26%). Land use linearly decreases from 1.02 ha per cow in the smallest farms to 0.68 ha in the largest farms (see Table 6). To produce 1 kg of raw milk the smallest farms require 1.71 m2 (including pastures), but the largest farms require only 0.76 m2 of land. The main contributors to the land use in the smallest farms are grass (50%) and protein feeds (27%), but as the farms intensify the role of grass decreases. Therefore, for the largest farms, the main contribution to land use is from protein pea (29%), grass (21%), alfalfa-grass silage (17.5%), and oat (14%). The method used to assess land use does not specify different land use categories, such as more or less productive agricultural lands.

Water footprint includes green (rainfall) and blue (freshwater stored in lakes, rivers, and aquifers) water use. The water footprint per kg of milk significantly decreases as farms get larger because the main contributor to water footprint is green water for grass production—which ranges from 88.6% in the smallest farms to 59.6% in the largest. Bluewater footprint (mostly used as drinking water for cows) accounts only for 1.7% of the water footprint in the smallest farms, increasing to 6.3% in the largest.

The differences in terrestrial acidification among the different size farms are minor (the coefficient of variation is 8%; see Table 4), but the impacts are highest for the farms with 10–50 cows and decreases for the largest farms. The largest contributor to terrestrial acidification for all the farms is ammonia emissions accounting for 53–67.5% of all impacts, but the next most important process is the growing of grass in pastures on the dairy farms due to manure from cows which is responsible for 17–24%. As ammonia accounts for most of the terrestrial acidification, the rankings are likely related to the amount and ruminal degradability of the protein as well as a level of milk production. As a result, it would seem that the amount of acidification could be constant across herd sizes if diets and cow genetic potential were equalized.

In the case of marine and freshwater eutrophication, the largest environmental impacts are on farms with 10–100 cows. The differences in these parameters among different management groups are large, with the coefficient of variation being 22% and 28%, respectively. Nitrogen flows into marine waters is estimated to be from 0.0037 kg Ne per kg of raw milk from the farms with 51–100 cows to 0.0018 kg Ne in the largest farms (see Table 4). The main contributor to marine eutrophication from milk production on the larger farms is alfalfa-grass silage feed (46% of all the impact from the farms with 51–100 cows to 30% in the largest farms), followed by grass grown at the farm (18–26%). However, for the smaller farms, the main processes contributing to marine eutrophication are grass and oat, which combined contribute more than 70% of all impacts. With regards to freshwater eutrophication, the two management groups which include the smallest farms (up to 50 cows) have the highest impact—0.00023 kg Pe/kg of FPCM. However, the main contribution to freshwater eutrophication from all the farms is from grass, followed by oat and alfalfa-grass silage.

Additionally, for freshwater ecotoxicity the largest environmental impact is from farms with 10–50 cows. The largest contributors are concentrates (oat ~30% and pea (24–31%)) (see Table 4). Alfalfa-grass silage is also responsible for a significant part of freshwater ecotoxicity—contributing 18% of all the terrestrial ecotoxicity impact from the smallest farms to 13% from the largest farms.

The results of our study demonstrate that the production of one kg of milk needs 37–62 g of oil equivalent fossil resources (see Table 4). In this case, the highest fossil resource intensity is for farms with 51–100 cows. The main contributors are production of alfalfa-grass silage (13–26%), protein feeds (17–25%), and fertilizer (P) (10–15%).

Additionally, we looked at changes in the environmental impacts of supplementary protein substitution. In this case, we analyzed substituting pea in the combined feed with soy, as it is increasingly entering the market as a high-protein feed for dairy cows in Latvia. For this, we used global process data on the environmental impacts of soy meal production from the Ecoinvent database [21]. Results demonstrated that most of the environmental impacts would increase, e.g., GWP for the farms with 10–50 cows would increase by 57% and for the largest farms by 32% (see Figure 2), but on average GWP per kg of milk would increase by 43%. As a result of the substitution, the land use in the smallest farms would increase by 12% and for other farms by 19%, on average increasing the land use of raw milk production by 18%. These increasing impacts are mostly linked to the soy productionrelated land-use change, domestic transport, and industrial processing [28]. However, it can also be assumed that introducing soybean meal to the diets of the dairy cows would decrease the protein supplement fed to keep the diets isonitrogenous.

**Figure 2.** Impact on land use and GWP per kg of FPCM if pea in cow diet was substituted by soy.

### **4. Discussion and Conclusions**

Our results show that the environmental impacts of milking cows are different on farms with different management strategies. For some of the impact categories, e.g., global warming, marine eutrophication, terrestrial acidification, and ecotoxicity, the highest environmental impacts are for mid-sized farms with 10 to 50 cows. In general, the trend of most environmental impacts is decreasing as the farms increase in size, thus partly confirming that many of the environmental impacts assessed here decrease with increasing milk yields [29–31]. However, the smallest farms perform worse than the mid-sized farms with 10–50 cows only in two out of eight environmental impact categories.

The largest farms (>200 cows), through increasing protein content in their feed and choosing more productive breeds, have managed to significantly increase their milk yields and thus, also improve the environmental performance on certain environmental impact categories of production. Our results demonstrate that dairy system intensification generally provides a significant decrease in all of the impact categories assessed, when measured per kg of raw milk. The most significant reduction was observed for water and land use. These results are somewhat different from others that argue that intensification decreases eutrophication and land use per kg of raw milk but has no clear effects on other impact categories.

In our study, for all the farms enteric fermentation had the largest impact on the GWP, followed by feeds (alfalfa, protein feeds, and grass) and electricity consumption. Increasing the consumption of concentrated feeds in the largest herds lead to higher milk yields and thus, helped to decrease relevant environmental impacts. The global warming potential of the raw milk at the farm gate in different LCA studies varies significantly—from 0.9 to 2.4 kg CO2e per kg of milk [27,32]. Results from our study are on the bottom line of this spectrum. This could be explained with the scope of the study excluding heifer rearing, which could increase the global warming potential of the milk by roughly 1/5 [33].

Differences in the environmental impacts of different farms can be explained by differences in feeds and diets of lactating cows as well as by differences in energy use of these farms, as the smallest farms are less energy-intensive in their performance but are more reliant on grazing (grass consumption). Higher milk yields are obtained in the larger farms, which have been achieved by incorporating proteins that are more resistant to ruminal degradation—soy or pea protein may be more efficiently used than proteins in forages like alfalfa.

It should be also noted that dairy cows in traditional systems provide other essential services, e.g., biodiversity protection and the conservation of cultural landscapes [34]. Unfortunately, as in the current study, these grazing benefits (e.g., carbon sequestration, biodiversity) are excluded due to the lack of approved assessment methods, and thus intensive systems are likely favored [35]. As a result, our study highlights the potential conflict between GHG emission and land use efficiencies and other environmental objectives, similar to Edwards-Jones et al. [36], and further emphasizes the need for the development of scientifically robust methods to assess carbon sequestration, biodiversity, and social impacts. Carbon sequestration is not covered by this study, but it is clear that optimized grazing, pasture improvement, and the restoration of degraded pastures are critical in increasing grassland carbon sequestration [37].

In addition, the specification and characterization of land and water use by different areas and crops could provide more insights. Grass production often takes place in areas unsuitable, or at least economically not profitable, for other crop production. Thus, the higher land use of small farms (of which 50% is of grass) may not reflect the higher land use of highly productive agricultural land. Similarly, the method used to assess water use in the paper does not reflect water availability or scarcity. Thus, also regarding water use, the higher impact of small farms may not reflect an actual higher water consumption of scarce water resources.

We can also assume that smaller farms tend to use cultivated local protein feeds which have a relatively lower environmental impact than the soybean byproducts commonly used in larger conventional farms. From this study, the total national pea consumption is estimated to be 157 thousand tons of peas a year. According to the FAO data in 2017 only 29.5 thousand tons of dry pea were produced in Latvia [38]. However, not all of the on-farm production of peas for domestic consumption is included in the national statistics.

It should be also considered that soy (one of the most popular high-protein feeds) production results in other social-ecological impacts [39], e.g., land acquisitions [40] and biodiversity loss [41], which are not properly covered by this study. In assessing soy production-related climate impacts, it would be also important to include LUC-related emissions. Additionally, local characterizations of water use to reflect a scarcity of water would improve the understanding of the water use of different crops from different countries.

One of the main policy recommendations from this study is to differentiate the policy approach to different size dairy farms in Latvia. The largest problems are with the management of farms with a herd size of 10–50 cows as they have some of the highest environmental impacts because their feeding strategy does not deliver higher milk yields. These farms should receive better training in feeding strategies as most likely they are overfeeding their cattle. Farms with a diverse range of cow breeds should also develop diversified feeding strategies for each of the breeds. The use of an automatic feeding system normally linked to milking robots could help to do this, as they increase the possibility of controlling the nutritional value of feed concentrates for each particular cow.

The smallest farms (≤9 cows) on another hand should be supported as self-sufficient providers and as important mechanisms for landscape and biodiversity conservation. Small farms with a higher fat content in the milk they produce should be encouraged to focus on differentiated production, e.g., cheese making. The smallest farms also should be supported in the use of cultivated local protein feeds and utilize cow manure effectively to create regenerative benefits in their field area [42]. Small farms are also important in terms of the efficiency of utilizing local sources of dry matter. Their feeding strategies appeared to be the most flexible of all the management groups. Thus, small farms also need support for this purpose. Diverse feeding strategies in small farms also would be a proxy for enhanced biodiversity.

However, we should also acknowledge that this study has several limitations linked to national data availability (e.g., national data on feed production, heifer breeding, differences among farms regards soil type, manure management, and the proximity to marine or aquatic habitats) and methodological shortcomings (e.g., excluding emissions of carbon

sequestration, the assessment method of water use, the use of proxy allocation, and the integration of social and biodiversity impacts in LCA). There is a need for the development of a national LCA inventory to further improve the quality and possibilities for environmental impact assessment.

**Author Contributions:** J.B. led the conceptualization, data collection, model development, calculations, interpretation of results, and writing of the manuscript; S.K. and H.H. contributed to research design, model development, analysis, interpretation of results, and writing of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partly financed by the specific support objective activity <IP_ADDRESS>. "Postdoctoral Research Aid" of the Republic of Latvia (Project No. <IP_ADDRESS>/VIAA/1/16/065 "Developing New Tools for the Sustainability Assessment of the Bioeconomy"), funded by the European Regional Development Fund (Project id. N. <IP_ADDRESS>/16/I/001). The APC was funded by authors.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

### *Article* **LCA to Estimate the Environmental Impact of Dairy Farms: A Case Study**

**Sara Zanni 1 , Mariana Roccaro 2,\* , Federica Bocedi 3 , Angelo Peli 2 and Alessandra Bonoli 4**

**Abstract:** Intensive farming is responsible for extreme environmental impacts under different aspects, among which global warming represents a major reason of concern. This is a quantitative problem linked to the farm size and a qualitative one, depending on farming methods and land management. The dairy sector is particularly relevant in terms of environmental impact, and new approaches to meeting sustainability goals at a global scale while meeting society's needs are necessary. The present study was carried out to assess the environmental impact of dairy cattle farms based on a life cycle assessment (LCA) model applied to a case study. These preliminary results show the possibility of identifying the most relevant impacts in terms of supplied products, such as animal feed and plastic packaging, accounting for 19% and 15% of impacts, respectively, and processes, in terms of energy and fuel consumption, accounting for 53% of impacts overall. In particular, the local consumption of fossil fuels for operations within the farm represents the most relevant item of impact, with a small margin for improvement. On the other hand, remarkable opportunities to reduce the impact can be outlined from the perspective of stronger partnerships with suppliers to promote the circularity of packaging and the sourcing of animal feed. Future studies may include the impact of drug administration and the analysis of social aspects of LCA.

**Keywords:** LCA; climate change; agro-livestock sector; GHG emissions; dairy farming

### **1. Introduction**

Climate change, described as the long-term heating of our planet caused by human activities since the pre-industrial period, is among the main challenges on a global scale [1]. As it is well known, one of the main causes of increasing temperatures is greenhouse gas (GHG) emissions in the atmosphere.

GHG concentration affects the global temperature by absorbing strongly radiant electromagnetic energy at wavelengths capable of emitting quantities of elevated heat [2].

According to their global warming potential, the main gases responsible for the greenhouse effect are methane, carbon dioxide, used as a reference for the phenomenon, water vapor, and nitrous oxide, each with different global warming potential and specific sources. Methane, for example, is about 21 times more efficient in trapping heat in the atmosphere compared to CO2, while nitrous oxide is 296 times more efficient, with persistence in the atmosphere for up to 114 years [3]).

The agro-livestock sector is responsible for a relevant share of these gases due to direct and indirect emissions, accounting for about 17% of the global GHG emissions in 2018 [4].

In addition to the release of GHG into the atmosphere, farming contributes in part to the worsening of air quality through the production of mainly nitrogen compounds from manure, particulates obtained from combustion, and volatile organic compounds

**Citation:** Zanni, S.; Roccaro, M.; Bocedi, F.; Peli, A.; Bonoli, A. LCA to Estimate the Environmental Impact of Dairy Farms: A Case Study. *Sustainability* **2022**, *14*, 6028. https:// doi.org/10.3390/su14106028

Academic Editor: Rajeev Bhat

Received: 5 April 2022 Accepted: 13 May 2022 Published: 16 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

other than methane (NMVOC), carbon black (BC), heavy metals (i.e., chromium, copper, nickel, selenium, zinc, lead, cadmium, mercury), dioxins, polychlorinated biphenyls (PCBs), polycyclic aromatic hydrocarbons (PAHs), hexachlorobenzene (HCB), and of particulates, both less than 10 microns (PM10) and 2.5 microns (PM2.5), deriving from several farm-level operations [5].

Livestock production has a significant environmental impact as it affects several natural resources, including land and soil, water, air, and biodiversity. Nevertheless, growing populations and economics have led to increased demand for animal products and the consequent expansion of the livestock sector over the past decades. At the same time, in developed countries, consumers demand animal products that are both animal welfare- and environment-friendly [3].

In 2019, world milk production grew to about 852 Mt and was forecast to grow at 1.6% per year over the next decade, faster than most other main agricultural products [6]. Globally, cattle are the largest contributors to total livestock greenhouse gas (GHG) emissions, producing about 65% (4.6 Gt CO2eq annually) of sector emissions, with milk production contributing 20% of total sector emissions [7]. Besides direct emissions, livestock systems are responsible for indirect emissions arising from land-use change, fertilizer use, energy, and transport emissions related to livestock operations and supply chains [8]. The three largest sources of GHG from milk production are emissions from manure management (CH4 and N2O), emissions from feed production, processing, and transport (CO2 and N2O), and emissions from enteric fermentation (CH4), the latter accounting for more than half the total of emissions [9].

Although absolute emissions from the dairy sector have increased in the last decades due to production growth in response to the increasing demand, dairy farming is becoming more efficient considering that emissions per unit of product are decreasing [8].

Strategies to reduce GHG emissions include changes in feeding, breeding, and management practices, which essentially lead to the intensification of livestock farming. Strategies that aim to increase productivity are very promising ways to reduce environmental impact; however, in most cases, they are likely to negatively impact animal welfare. For example, intensive housing conditions increase the risk of social stress or hinder the expression of natural behavior [10].

The rising GHG emissions require shifting production systems toward carbon neutrality in order to take action to combat climate change and its impacts, which is one of the United Nations 2030 agenda for sustainable development goals (SDG 13). Moreover, sustainable food production systems and resilient agricultural practices are among the targets to achieve the "zero hunger" goal (SDG 2) and to promote new consumption and production models (SDG 12) by changing the way we produce and consume, mainly in the agri-food sector where these elements are so close, in the farm-to-fork perspective. At the same time, progress toward the sustainable use of natural resources is key to protecting biodiversity and ecosystems (SDG 15). Improving animal health and welfare and reducing GHG emissions through new farming techniques to meet societal demands on safe and sustainable food and contribute to climate change mitigation are among the nine key objectives of the 2023–2027 European common agricultural policy.

The relationship between livestock production and climate change is two-sided. Climate change, deriving from global warming, adversely impacts livestock production, with direct effects on animal health, reproductive efficiency, production performances, and behavior, but also indirect effects deriving from changes in the quality and quantity of feed, water availability, ecosystem alterations leading to changes in the biology and distribution of pathogens and vector-borne diseases [11]. An increase in temperatures negatively affects milk production, especially in cows with higher milk yield and milk quality, with a reduction in casein content [12].

Dairy cattle production systems need to adapt to climate change, but, on the other hand, they must commit to contributing to GHG reduction targets and minimize other negative environmental impacts while continuing to meet society's needs. Agriculture is estimated to generate 11% of all global emissions [13], and dairy farm contribution to the overall impact of milk products is almost as high as 72% [14]. A special focus on this production phase is required to meet the global goals.

The dairy sector is an extremely complex system with numerous interacting components; consequently, determining the best strategies to reduce GHG emissions will depend on each farm's local conditions and objectives. Therefore, a deeper understanding of the underlying driving factors of dairy cattle farming environmental impact is highly relevant to support strategies for limiting adverse effects on the environment and protecting the livestock sector [15].

Studies on the environmental impact assessment of dairy products have increased dramatically in the last decade and have been reviewed by several authors [16–21]. The main goals of these studies were either to assess the potential environmental impact of the product [22–25] or to compare different management systems [26–28]. However, very few of these have focused on Italian production systems [29–32] and even less on Italian Protected Designation of Origin (PDO) dairy products (e.g., Parmigiano Reggiano, Grana Padano) [33]. In Italy, milk production is concentrated in the northern plain regions, characterized by little rain during summer and a highly urbanized territory with low availability of agricultural land; cattle are kept mostly indoors, and maize silage is the main forage crop. Moreover, PDO product specifications impose restrictions on management strategies, sourcing, and type of animal feed.

Considering the peculiarity of the environmental conditions and farming practices, this study aimed to set up an LCA model specific for intensive dairy cattle farms involved in PDO production in Northern Italy. We aimed to evaluate the most relevant items of environmental impact triggered by this activity, considering its ability to impair animal health and welfare. To this end, a preliminary assessment was performed on a dairy cattle farm in the Emilia-Romagna region in Italy. The LCA model was realized using Simapro 8.5.1, including farm processes and main supplies, in terms of animal feed and end-of-life processes of the most relevant waste flows.

With the ultimate aim of integrating farm environmental impact and animal health and welfare, an interdisciplinary research approach combining engineering and veterinary expertise was carried out to evaluate farms operating conditions in terms of animal management and welfare, and specifically, to assess environmental impacts based on an LCA approach [34]. This approach offers grounds for future research and development of targeted sustainability strategies in the sector, specific to the studied geographical area, with a concrete possibility of application in the field.

### **2. Materials and Methods**

The following sections report a detailed description of the site of the study (Section 2.1) and the conceptual model developed for the implementation of LCA (Section 2.2).

### *2.1. Site Description*

The present study was carried out on an Italian dairy cattle farm producing Parmigiano Reggiano and soft cheeses in the Emilia Romagna Region.

Prior to the farm visit, a form for primary data collection was designed. During the farm visit (July 2020), data were collected on paper and subsequently transferred to Microsoft® Excel files. Data sources included the farm management software records, information collected by interviewing the farmer, and invoices.

For data that could not be collected on-farm, default values were defined based on processes available on Ecoinvent v.3.5 [35], with geographical reference to the European context, where possible, or Italian specific data, as in the case of the electricity production mix.

Information on farm size, number of animals, productivity, animal welfare, drug use, feed management, water consumption, fuel consumption, electric energy, manure management, and plastic waste were collected (as summarized in Table 1). It was not possible to quantify the use of chemical fertilizers and bedding material.

**Table 1.** Farm-specific primary data used to design the conceptual model for the assessment.

The farm is structured in four sites, located between Bologna and Modena provinces. Site A hosted calves and replacements up to 5 months of age, pregnant heifers (from two months pregnancy to delivery), and dry cows. Replacements from 5 months of age to the first month of pregnancy were hosted in site B. Site C hosted mostly primiparous milking cows, and site D hosted mostly multiparous milking cows. Animals are moved from site A to site C or D twice a week and from site A to site B twice a month.

The most represented breed was Holstein Friesian, but Italian Red Pied, Jersey, and Montbéliarde x Swedish Red and White crossbreds were also present.

Average daily milk production amounts to 35.15 kg for multiparous cows and 30.15 kg for primiparous cows. The milk is destined to produce Parmigiano Reggiano.

Cows are loose housed on deep litter with straw bedding (dry cows and heifers) or in pens with cubicles covered with sand (milking cows). Dry cows and heifers have access to external exercise areas. The facilities are equipped with ventilation systems consisting of ceiling fans for the resting areas and horizontal flow fans for the feeding zones. The waiting parlor is equipped with a ceiling fan and water sprayers.

Calves up to one month of age are reared in individual pens and then moved to group pens. They are fed through an automatic calf feeder and weaned at 70 days of age.

Site C and site D scored 70.33/100 and 71.31/100, respectively, for animal welfare in 2018, as certified by the Italian National Reference Center for Animal Welfare (CReNBA).

The diet consists of 60% fodder and 40% concentrates. Fodder (alfa-alfa, grass hay) and part of sorghum are self-produced; the other feeds are bought in. Given the production requirements for Parmigiano Reggiano, corn silage is only fed to heifers (site B). Daily average dry matter intake is 25 kg for milking cows, 15 kg for dry cows and pregnant heifers, and 8 kg for replacements.

The water comes from the aqueduct. Site C was also supplied from a well, but it was not possible to quantify the water used. There is no water recycling system in place.

The slurry is spread on the farm's land intended for fodder production (470 ha). Figure 1 reports the data collection at the basis of the conceptual model definition. *Sustainability* **2022**, *14*, x FOR PEER REVIEW 6 of 15

**Figure 1.** Conceptual model flowchart. **Figure 1.** Conceptual model flowchart.

#### The LCA methodology was applied based on ISO 14040 [34] using the Simapro 8.5.1 *2.2. Conceptual Model for the Assessment*

software and the Ecoinvent 3.4. database [35]. The Functional Unit (FU) was set as 1 kg of Fat and Protein Corrected Milk (FPCM) produced in 2019, in accordance with similar studies in the field [24,25,33], and it is defined by [36] as FPCM (kg) = raw milk (kg) × (0.337 + 0.116 × fat content (%) + 0.06 × protein The conceptual model at the basis of the LCA has been defined based on a cradle to farm gate perspective over a calendar year. Considering the specificity of the farm analyses, the model has been developed based on primary data, whenever possible, and previous studies in the field (e.g., Grassauer et al., 2022 [28]).

content (%)) The average Italian national fat and protein contents reported in the year prior to the Several assumptions have been introduced aimed at reducing the uncertainties of the results and at providing a preliminary assessment of the generated impacts, as presented below:

**3. Results and Discussion**  In developed countries, since more and more citizens are concerned about health, A specific focus was dedicated to animal feed, which has been represented in the following terms, considering each upstream process and relative impacts:

environmental, ethical, and animal welfare issues, consumption of animal products has

2050 and 11.2 billion by 2100 [40], while fulfilling the obligation to reduce greenhouse gas

Despite the negative effects on the environment, our planet faces the challenge of

For these reasons, the involvement of milk and dairy producers at the local level and