INTENSIFICATION OF GRASSLAND-BASED DAIRY PRODUCTION AND ITS IMPACTS ON LAND, NITROGEN AND PHOSPHORUS USE EFFICIENCIES

Many grassland-based dairy farms are intensifying production, i.e., produce more milk per ha of land in response to the increasing demand for milk (by about 2% per year) in a globalized market. However, intensive dairy farming has been implicated for its resources use, ammonia and greenhouse gas emissions, and eutrophication impacts. This paper addresses the question of how the intensity of dairy production relates to N and P surpluses and use efficiencies on farms subjected to agri-environmental regulations. Detailed monitoring data were analyzed from 2858 grassland-based dairy farms in The Netherlands for the year 2015. The farms produced on average 925 Mg·yr − 1 milk. Milk production per ha ranged from < 10 to > 30 Mg·ha − 1 ·yr − 1 . Purchased feed and manure export strongly increased with the level of intensification. Surpluses of N and P at farm level remained constant and ammonia emissions per kg milk decreased with the level of intensification. In conclusion, N and P surpluses did not differ much among dairy farms greatly differing in intensity due to legal N and P application limits and obligatory export of manure surpluses to other farms. Further, N and P use efficiencies also did not differ among dairy farms differing in intensity provided the externalization of feed production was accounted for. This paper provides lessons for proper monitoring and control of N and P cycling in dairy farming.

Please note that technical editing may introduce minor changes to the manuscript text and/or graphics which may affect the content, and all legal disclaimers that apply to the journal pertain. In no event shall HEP be held responsible for errors or consequences arising from the use of any information contained in these "Just Accepted" manuscripts. To cite this manuscript please use its Digital Object Identifier (DOI(r)), which is identical for all formats of publication.

INTRODUCTION
Global consumption of dairy products is projected to increase by about 2% per year and mean consumption per capita by about 1% per year during the period 2018-2027 [1] . Most of the increase will take place in Asia and Africa. Mean per capita milk consumption was about 110 kg during 2015-2017, but with large differences among countries. Total production was almost 850 Tg in 2017; about 83% was cow milk, 12% buffalo milk and 5% was from sheep, goat and camels (IDF, 2018). Most of the production was in Asia (30%), followed by the European Union (EU-28; 28%), North and Central America (18%) and South America (9%). The self-sufficiency was 90% for Asia, 109% for EU-28, 100% for North America, 84% for Africa and 290% for Oceania in 2017 [2] .
The projected 25% increase in global consumption of dairy products (fresh milk, yogurt, milk powder, butter and cheese) between 2018 and 2027 [1] raises the question of how and where to produce the additional milk [3] . This question is not only relevant for Asia and Africa, where most of the consumption increase will take place and where the self-sufficiency is < 100%, but also for the main dairy exporters Oceania, EU-28 and North America, where dairy farming is under pressure to reduce greenhouse gas and ammonia emissions, and N and P losses [4][5][6] . This suggests that producer price is not the only factor in the global market that determines where the additional production will take place, despite the fact that the main dairy exporters currently have the lowest producer prices [2,7] . The aforementioned question has to be considered also in the context of increasing consumption (and hence production) of other food commodities (cereals, oil, meat, fish, vegetables and fruit) during the next decade (by 1%-2% per year), while the total area of arable land and pasture will not increase [1,8] . As a consequence, resource use efficiency in agriculture will have to increase, including in dairy production [9][10][11] .
Intensification of dairy production is commonly expressed as increasing milk yield per ha of land, per farm or per laborer. Intensification of dairy production is a complex process, as it has many driving forces, biophysical, socioeconomic and environmental constraints, and often unintended consequences. In general, intensification is a result of technological progress, which is driven by developments in markets, technology and/or policy. Sustainable intensification is commonly defined as 'an economically profitable, socially acceptable and environmentally sound production increase, in balance with a responsible consumption, for now and later' [12,13] . Sustainable intensification is commonly seen as the way forward, although there is debate about the actual meaning of the various words in the aforementioned definition, about the balance between economic, social and environmental requirements, about responsible consumption and the balances between here and there, and now and later [13] . Sustainable intensification in dairy production commonly encompasses the notion that the increase in milk yield must be achieved with less inputs and with less pollution of the environment, while protecting birdlife, biodiversity, and natural landscape elements [14,15] .
There are millions of dairy farms in the world encompassing a bewildering diversity of systems in a wide variety of socio-economic and environmental conditions [2,16] . This diversity and complexity are often not considered fully in global scoping and assessment studies. Economically competitive farmers are often well-educated and have easy access to markets, capital, technology and advice [17] , which makes farmers in affluent countries often quite prosperous [18] . These farmers are, however, not necessarily better managers of resources, including nutrients. Intensification of dairy production, i.e., increasing milk yield per cow, per ha of land and/or per laborer, and upscaling, i.e., more land and cows per farm, are common strategies for dairy farms to increase competitiveness and to stay in business [19,20] , although some farmers diversify and seek additional income sources [21] . Intensification and upscaling at farm level are commonly associated with changes in production structure, resource use efficiency and emissions per ha of land and per kg of milk produced. It may lead to outsourcing of certain activities, including animal feed production and manure disposal at other farms. The relationships between intensification and upscaling of dairy farms on the one hand and resource use efficiency and emissions per ha of land and per kg of milk produced on the other hand are quantitatively not well understood. This lack of understanding hinders extension services and policymakers in guiding dairy farms adequately in the development process.
This study aimed at increasing the quantitative understanding of the relationships between the intensity of milk production and resource use efficiency and emissions per ha of land and per kg of milk produced level on grassland-based dairy farms subjected to agri-environmental regulations. We analyzed data from dairy farms in The Netherlands, which are among the most intensive and productive grassland-based dairy farms in Europe [22] and the world [23] . There were 17,000 dairy farms which produced 14.1 Tg of milk in 2018 [7] . These farms are required to adopt agri-environmental regulations depending in part on the intensity of milk production [24] . Dairy farms used about two-third of the total agricultural area (1.8 Mha) in The Netherlands, and in addition import feed from other countries. We statistically analyzed data from > 5000 farms for the years 2013-2015 but focus on data from 2858 farms for the year 2015.

Dairy farming characteristics in the Netherlands
Most dairy farms are family farms. In 2018, these farms had on average 95 dairy cows and 64 ha of land (note, median numbers are smaller than average numbers). Farms with high milk production per ha of farmland were mostly situated on sandy soil where farm size is relatively small compared to the farms on clay soils for historic reasons (households had more children and farms were smaller on 'poor' sandy soils than on 'fertile' marine and riverine clay soils [25] ). On average 80% of the farmland is managed grassland and the other 20% is used for (silage) maize production. Many dairy farms exchange part of the land for one year with arable farmers to grow potatoes or bulb flowers (mainly tulips and lilies), to widen the crop rotation of arable farms, and receive 2000-3000 EUR·ha −1 ·yr −1 in exchange. Purchase costs of land range from 40,000 to 100,000 EUR·ha −1 , depending on the region and the land quality. Larger farms (> 150 cows) may have one or more laborers and most dairy farms make use of professional contractors for field work (manure application and harvesting). About 25% of the dairy farms make use of automatic milking systems [7] .
The average ration of the dairy herd including young stock consisted of 20% fresh grass, 35% grass silage, 15% whole maize silage and 30% concentrates (in dry matter). About 80% of the farms used rotational grazing during the summer half year from May to October, often during the daytime only (mean grazing duration was 1650 h·cow −1 in 2018). Mean milk yield was 8800 kg·cow −1 ·yr −1 and 13,000 kg·ha −1 in 2018. Before 2015, total milk production was limited by tradable milk quota, and from 2017 by tradable production rights, expressed in kg P2O5 (1 kg P2O5 is equivalent on average to 200-230 kg milk, depending on milk production efficiency). The purchase costs of the production rights depend on the market, but were in the range 200 to 250 EUR (kg P2O5) −1 in 2019 (roughly 8000-10,000 EUR·cow −1 ) [7] .
Manure and fertilizer use are constrained by N and P application limits [24] . In 2015, manure N application limits ranged from 170 to 250 kg·ha −1 ·yr −1 , depending on farm-specific permits. Total P application limits ranged from 17 to 52 kg·ha −1 ·yr −1 for arable land and from 33 to 52 kg·ha −1 ·yr −1 P for grassland, depending on soil P status. In practice, P application limits define the maximum manure application rates on farmland with a high soil P status, while manure N application limits define the maximum manure application on farmland with low soil P status. By default, all farmland is assumed to have a high soil P status, unless farmers can show through soil P tests by accredited laboratories that the soil P status of the farmland is low to moderately high (there are 5 classes). There is a ban on the use of synthetic P fertilizer on most dairy farms [24] .
Most of the dairy manure (95%) of housed cattle is collected as slurry below slatted floor in cubicle houses. A few farms have littered barns and solid manure. The slurry has to be stored in leakproof and covered manure stores, and has to be applied to land using low-emission techniques during the growing season only. Manure production, manure use on farmer's own land, and manure export have to be registered at a governmental office, which verifies the submitted data. Farms with a surplus of manure have to process and export the surplus manure to other farms, at a cost of 10-20 EUR·m −3 of slurry (one dairy cow produces 25-30 m 3 ·yr −1 ). Slurry separation allows for cost-effective export of a relatively P-rich and dry solid fraction in a small volume at reduced cost. The N-rich liquid fraction may then be landapplied on the farm's cropland within the set manure N application limits [24] .

Data collection
Data from dairy farms were collected in the test phase (2013-2015) of the KringloopWijzer model, which is a nutrient management and assessment tool, and a joint initiative of Wageningen University & Research, the dairy sector, feed suppliers, processing industry, farm advisors and accountants, and the Ministry of Agriculture, Nature and Food Quality. Dairy farms participated on a voluntary basis in the test phase. Farmers were guided by advisors to complete the data recording in the KringloopWijzer model, and use the results for management decisions. All (input) data and results are stored in a central online database, owned and supported by the dairy industry. Permission was obtained from the participating organizations to analyze the data from the test phase.
The origin of the KringloopWijzer model is based on long-term research at experimental dairy farm De Marke [26,27] and on monitoring data of pilot commercial dairy farms [28] . The model has served as a learning and management tool, also through exchange of information among farmers in study groups [29] . The model has been extensively tested through experimental measurements and uncertainty analyses [30,31] , before large-scale use in practice.
For the assessment of N and P flows and use efficiencies of dairy farms, four components are distinguished within the dairy farm, i.e., herd, manure, soil and crop [32,33] . The KringloopWijzer model estimates the N and P flows through these components, i.e., an output from one component is an input into another component, but losses are incurred in these transfers. The N and P balance of a component, i.e., the difference between inputs and outputs, characterizes the (in)efficiency in management of a particular nutrient in a particular part of the farm, revealing the weakest and strongest parts of the farming system [34] . Hence, N and P balances were assessed at farm level and at component level.
The KringloopWijzer model and the data collection procedures are described in detail in the Supplementary Information ( Fig. S1 and Table S1). All annual N and P inputs and outputs at farm level were recorded, based on farm accounts (purchased feed, fertilizers, manure, and exported milk, animals, crop products and manure) or literature data (atmospheric N deposition and biological N2 fixation). Total annual N and P surpluses at farm level were derived from the difference in total input and total outputs. The N and P use efficiencies (NUE and PUE) were derived from the ratio of total output and total input. Note that manure export was recorded as a negative import (and not as an output), following Quemada et al. [22] , because this way of recording provides estimations of NUE and PUE which allow a more precise comparison among dairy farms greatly differing in intensity and manure export.
The farm-internal NUE and PUE of the herd ( Fig. S1 and Table S1) were based on the ratio of the N and P outputs in milk and meat and the N and P inputs in feed. The NUE and PUE of the soil component was based on the ratio of the N and P outputs in harvested crop (including through grazing) and the N and P inputs in manure, fertilizers, atmospheric N deposition and biological N2 fixation. The yields of homegrown maize silage and herbage (both fresh grass and silage) were estimated by difference; i.e., the energy requirements of the dairy herd (as function of breed, age, milk production, lactation, pregnancy, and housing/grazing) minus the sum of the energy in imported feed. Corrections were made for changes in the stocks of maize silage and grass silage. The relative proportions of the intake of fresh grass during grazing, silage grass, maize silage and concentrates by the dairy cattle were based on estimations of the available feed stuff on the farm by the farmer, but checked and eventually revised following calculations of the feed energy balance of the whole farm, using an optimization routine. Samples of grass silage and maize silage were analyzed on most farms, by accredited laboratories on farmers' request. The estimated homegrown yields of maize silage and herbage have been tested on a sample of farms during a period of 10 years and the accuracy confirmed to be within 5% (see Supplementary Information, Fig. S2). A detailed description of the algorithms of the KringloopWijzer model can be found in De Vries et al. [35] Table S2 provides a summary overview of the farm data from the test phase (2013-2015). All data were checked for plausibility, using a checklist with criteria. Following this initial data screening, a sample remained of 1096 farms from 2013, 1597 farms from 2014 and 2858 farms from 2015. Main reasons for rejection of farm data (about 15% of the farms) were (1) incomplete records, and (2) unlikely records (e.g., harvested grass yields of > 20 t·ha −1 DM). We estimated a wide range of indicator values and found that the mean indicator values were rather similar for the three years. Based on the similarity, we decided to analyze the data from 2015 in more detail because the number of farms from 2015 was larger than those from the earlier years, and because of the learning effect; many farms in the sample from 2015 also participated in the sample from 2013 and 2014, and fewer farm data sets had to be excluded from the analyses due to questionable input data. The sample of farms is not meant to be representative for dairy farms in The Netherlands; the average farm size and intensification level was slightly higher than the mean of all dairy farms. Yet, the wide range of volunteering farmers in the test phase does all for examination of the relationships between the intensification level of dairy production and NUE and PUE, and related impacts.

Data analysis
Farms were ordered in ascending milk production per ha, and then divided in 10 groups of equal size (286 farms per group). The variation in indicator value for each group of farms was visualized by means of box plots, and statistically analyzed by means of ANOVA to test whether there are differences between groups. ANOVA was followed by pairwise testing, employing a significance level of 5%.
Following Quemada et al. [22] , the effect of externalization of the N and P cost of feed production, through the purchase of feed from elsewhere, was explored by assuming that the purchased feed was produced with a NUE and PUE of 50% and 75%. Hence, we increased the N and P inputs from purchased feed by a factor of 2 and 1.33 in additional sensitivity analyses.
Relationships between explanatory indicators and NUE at farm level were analyzed with multiple regression models and correlations (Tables S3 and S5). We established regression models for NUE_farm using a range of potential explanatory indicators (regressors). The relevance of variables was tested using the RSEARCH procedure in GenStat (all possible subset selection) [36] . Only explanatory indicators (regressors) that were sufficiently uncorrelated (r < 0.70) have been included in the selection process to avoid the problem of collinearity [37] . In case of high correlations, one of the variables was selected for inclusion in the selection process and the other was rejected. To identify the best parameter combinations, the percentage of variance accounted for (R 2 adj; i.e., adjusted for the number of parameters), the value of Mallows' Cp [37] , and the p value of the parameter estimates were evaluated. We selected models with the highest R 2 adj and lowest Cp, and significant parameters (p < 0.05) (Table S4).

Intensity levels of dairy production
The distribution of mean milk production per ha of farmland was skewed to the right (Fig. S3). The variation was large; farms of group 1 had a mean milk production of 11 Mg·ha −1 ·yr −1 and group 10 farms had on average 30 Mg·ha −1 ·yr −1 (Table 1; Fig. 1). Farms with high milk production per ha of farmland also had a relatively high milk production per cow, a high feed protein import, but a relatively small area of farmland (Fig. 1). This indicates that farms had increased milk production per ha of farmland, if the area of farmland was relatively small. Intensive dairy farms had a low number of young stock (calves and heifers), and a low grazing intensity (Table 1), suggesting that they aimed at maximizing milk production per ha of farm land. Mean total milk production per farm ranged from 651 Mg·yr −1 for group 1 to 1191 Mg·yr −1 for group 10. Mean feed conversion efficiency increased significantly with intensification level; it ranged from 0.98 kg milk per kg feed intake for group 1 to 1.15 kg·kg −1 for group 10 ( Table 1). The increased feed conversion efficiency reflects the increased milk production per cow and the decrease in the number of young stock per cow with the intensity of milk production.
The variation among farms in feed protein import was notably large for group 10 (Table 1), which was related in part to the large variation in milk production per ha of farmland. However, there was also a significant variation within the other group of farms in feed protein import, which may be related to variations in milk production per cow and in land productivity. Mean dry matter yield of grass ranged from 9.7 Mg·ha −1 ·yr −1 for group 1 to 11.6 Mg·ha −1 ·yr −1 for group 10. Mean dry matter yield of whole maize plants ranged from 18.3 Mg·ha −1 ·yr −1 for group 1 to 17.7 Mg·ha −1 ·yr −1 for group 10 (Table 2). Variations in yield were in part related to soil type (not shown).
Dairy farms on sandy soils prefer to have a relatively large area of land for silage maize production because sandy soils are drought sensitive and maize (a C4 species) has a higher water use efficiency and a higher yield potential than perennial ryegrass (Lolium perenne) (a C3 species), which is the dominant grass in grassland in the Netherlands. However, governmental regulations limit the area of maize at farm level to 20% because of its nitrate leaching risks [38,39] . Hence, the shares of the area of grassland and maize land were remarkably constant across the 10 groups (Table 1). However, the share of silage maize in the diet of the herd of most intensive farms was higher than that of extensive farms because the shortage of feed at intensive farms was covered in part by the purchase of maize silage. Notably, the crude protein content of the diet of the herd remained largely constant across intensification levels, at 158 g·kg −1 (Table 2, Fig. S4).  Note: Farms were ordered in ascending order of milk production and then divided in 10 groups of 286 dairy farms each; means are presented per group. Values followed by the same letter within rows are not significantly different between groups (p < 0.05). Note: Farms were ordered in ascending order of milk production and then divided in 10 groups of 286 dairy farms each; means are presented per group. Values followed by the same letter within rows are not significantly different between groups (p < 0.05).

Nitrogen balances and use efficiencies
The N surplus at farm level increased slightly with intensification level (Fig. 2). Mean N surplus was 174 kg·ha −1 ·yr −1 for group 1 and 208 kg·ha −1 ·yr −1 for group 10 (Table 1). Linear regression analysis indicated that the mean N surplus increased by 1.5 kg·ha −1 ·yr −1 when milk production increased by 1000 kg·ha −1 ·yr −1 (Table S3). The relative constancy of the N surplus at farm level indicates that the increasing N inputs from animal feed (Fig. 1) were balanced by increasing N outputs in the sale of milk and animals (Fig. 2(c)) and the export of animal manure (Fig. 2(d)). The NUE at farm level strongly increased with intensification level, from a mean of 34% for group 1 to 49% for group 10. This strong increase in NUE is related to a number of factors, including the increase in milk production per cow, the decrease in the number of young stock per cow, the decrease in grazing intensity, the export of animal manure and the offsite (externalization of) feed production (see section 3.4), with increasing intensity.  (Table S4). Milk production per ha was strongly correlated with the purchase of feed, the sales of milk and animals, and the export of manure, and hence could not be combined in one regression model. A multiple regression model with (1) purchase of feed, (2) NUE of the soil/crop component, (3) NUE of the herd, and (4) the sales of milk and animals explained 83% of the variance (Table S4). The best model explained 84%; hence, additional factors did not contribute much.
The NUE of the herd increased significantly from a mean of 23% for group 1 to a mean of 26% for group 10 ( Table 2). The NUE of the soil/crop component increased from a mean of 66% for group 1 to a mean of 71% for group 10 (Table 2). These results are consistent with the increased feed conversion ratio mentioned above, and with the pressure of intensive farms to get as much feed from their own land as possible, respectively. However, there was a large variation within groups ( Fig. 2 and Fig. 3).
Manure N export strongly increased with milk production per ha (Fig. 2(d)). Farms have to export manure N and/or P obligatory, when total manure N and/or P production at farm level (corrected for gaseous N losses during storage) exceeds limits placed on application for manure N and/or P on the land of the farm (see section 2.1).

Phosphorus balances and use efficiencies
The P surplus at farm level decreased slightly with intensification level (Fig. 3). Mean P surplus was 1 kg·ha −1 ·yr −1 for group 1 and −7 kg·ha −1 ·yr −1 for group 10. Hence, average P surplus decreased by about 0.4 kg·ha −1 ·yr −1 when milk production increased by 1000 kg·ha −1 ·yr −1 , and the P surplus turned into a P deficit at a high level of milk production per ha. The relative constancy of the P surplus at farm level indicates that the increasing P inputs from animal feed (Fig. 3) were more than balanced by increasing P outputs in the sale of milk and animals (Fig. 3(c)) and the export of animal manure (Fig. 3(d)). Intensive dairy farms have to export relatively large amounts of manure P, also because many intensive farms have farmland with high soil P status, and thus low P application limits. Most of these farms have been intensive for some time and have enriched the soil with manure P in the past, when the N and P application limits were less strict [24] . Government regulations now require farmers to reverse the situation; P application limits for land with a high soil P status are set at such a level that annually 5-10 kg of P will be mined from the soil, until the soil P status drops to a level with higher P application limits. Our results indicate that soil P was mined at about 70% of the dairy farms (Table 1). Fertilizer P applications were on average ≤ 0.2 kg·ha −1 ·yr −1 , which reflects the ban on mineral P fertilizer application on most dairy farms.
The PUE at farm level increased with intensification level, from a mean of 120% for group 1 to 140% for group 10. These high PUE values reflect that total P output was larger than total P input, on basically all farms. Mean PUE was also higher than 100% in 2014, but not in 2013 (Table S2). Such annual variations likely reflect annual changes in feed stock and manure storage. Manure P export strongly increased with milk production per ha ( Fig. 3(d)). Manure P export increased from a mean of 2 kg·ha −1 ·yr −1 for group 1 to a mean of 38 kg·ha −1 ·yr −1 for group 10. In group 10, similar quantities of P left the farm in milk and meat as in exported manure (Table 1).
The PUE of the herd increased from a mean of 30% for group 1 to a mean of 34% for group 10 ( Table 2). PUE of the soil/crop component of the farm increased from a mean of 99% for group 1 to a mean of 122% for group 10 (Table 2). A PUE of the soil/cop component exceeding 100% reflects soil P mining.

Correcting NUE and PUE at farm level for externalization of feed production
When the imported animal feed N and P were corrected for the likely NUE and PUE of feed production (i.e., 50%-75%), the NUE and PUE at farm level were lower and increased much less with intensification level (Fig. 4). This suggests that the strong increases in NUE and PUE at farm level with the level of intensification ( Fig. 2(b) and Fig. 3(b)) is in large part related to offsite (externalization of) feed production. Assigning a value of 50% for NUE to purchased feed stabilized farm NUE at a level of about 22% across all 10 groups (Fig. 4(a)). Assuming a NUE of 75% for purchased feed still gave a significant increase in farm NUE, from 26% for group 1 to 35% for group 10 (Fig. 4(c)). In contrast, a value of 50% for PUE of purchased feed led to a significant decrease in farm PUE from about 45% for group 1 to 40% for group 10 ( Fig. 4(b)), while assuming a PUE of 75% for purchased feed stabilized farm PUE across groups at 75-80% (Fig. 4(d)). In summary, assigning an externalization production efficiency to purchased feed has large effects on farm NUE and farm PUE.
Assuming an efficiency of 50% of the imported feed Assuming an efficiency of 75% of the imported feed

Intensification of dairy production in the Netherlands
Dairy production and consumption have a long history in The Netherlands. The first domesticated cattle probably appeared some 7000 years ago in Western Europe [40] . Initially, cattle were used for meat production, but soon also for traction, leather and milk. Dairy cattle were kept on small mixed farms which also had pigs, laying hens and arable land for potato and cereal production, until the 1950s [41] . Governments started to support agricultural modernization from the end of the nineteenth century through stimulating education, research and extension, and later on also through structural measures, including land reclamation [25] and subsidies for farm buildings and animal breeding programs, although the latter were soon privatized [42] . Government support increased in the second half of the twentieth century through the establishment of the EU common agricultural policy in 1962. Also, markets became larger and improved technology became available. As a result, dairy production became specialized on grassland-based farms, milk production per cow and ha of land increased rapidly, and farms became larger while the number of farms decreased rapidly [43] .
The incentives for intensification of dairy production ultimately led to a milk surplus and to nutrient surpluses. In response to the milk surplus, a milk quota system was implemented in 1984, which regulated total milk production at national level and farm level, and indirectly stabilized the milk price. In response to the N and P surpluses, regulations related to manure and fertilizer management were implemented from 1990s onwards, which decreased the degrees of freedom in dairy farming [24,44] . The abandoning of the milk quota system in 2015, following agreements within the framework of the World Trade Organization, provided an incentive to intensify dairy production again, but also put increasing pressure on the manure management regulations [45][46][47] . In response, dairy manure P production rights were implemented from 2017.
Our farm data were from 2015, which is at the start of latest wave of dairy production intensification. The intensification in the short period in between the regulation by milk quota and the regulation by manure P production rights has contributed to the skewness of the distribution of milk production per ha farmland (Fig. S3). Further intensification is still possible, provided farmers are able to purchase the manure P production rights, have permits for animal houses, and house the additional cattle in low-emission barns, and export the surplus manure. In addition, farms with a manure surplus can only intensify further if additional land is purchased; the permissible (legal) manure application on this additional land must cover at least 50% of the additional manure production. In summary, further intensification is strongly constrained, both legally and economically. Only the most competitive farmers are able to intensify further, with the help of financial advisors because of the complexities involved.
About two-thirds of the total milk production in The Netherlands are exported, following processing, which occurs for about 80% by the farm cooperative FrieslandCampina. The processing of the milk and the export of milk specialties have a positive effect on the milk price, and have allowed the dairy sector to keep its position in the world market [2,7] . The influence of processing industries, retail and NGOs on dairy farming has increased during recent decades. For example, dairy farmers get a higher milk price if the performance of dairy production meets criteria related to the On the Way to PlanetProof concept [48] . There is also pressure from some NGOs, farmer organizations, and political parties to reward dairy farmers, through a higher milk price, when farms produce at least 65% of the feed need on their own farm, and contribute to a more circular economy [49] . These movements may limit further intensification of dairy production and make dairy production more farmland-based.

Impacts of the intensification of dairy production
A main driver for intensification of dairy production is the need for farmers to reduce cost and remain competitive in a global market [14,50,51] . The incentive to intensify has been greater on dairy farms with a relatively small farmland area; i.e., the most intensive farms had the smallest land areas (Fig. 1). Dairy farms with a large milk production per ha of farm land also had a relatively high milk production per cow ( Fig. 1(b)) and required large inputs of animal feed (Table 1). Fertilizer N input increased only slightly with intensification level and not further beyond a milk production of about 15.5 Mg·ha −1 ·yr −1 . Feed protein-N import increased on average by 18.5 kg per 1000 kg of milk (Table S3), which roughly translates to a mean feed-N into milk-N conversion of about 30% (using a mean milk protein-N content of 5.6 g·kg −1 ) of the purchased feed. An NUE for purchased feed of 30% is higher than the overall feed-N use of the herd (23%-26%; Table 2), likely because of differences in feed quality and protein content between purchased and farm-grown feeds, and because of the improved utilization of farm-grown feed with an increase in intensification level.
The NUE at farm level was positively related to the intensity of milk production; mean NUE was 34% for group 1 and 49% for group 10. This is in contrast with reports indicating that NUE of grassland-based dairy farms decrease and N surpluses increase with an increase in milk production per ha [52,53] . The decrease in NUE and increase in N surplus in these latter studies were mainly a result of the increased N fertilizer input (to increase homegrown feed production), while there was no manure export at a high level of milk production per ha. Further, the NUE at farm level was relatively high compared to NUE values of dairy farms in other countries [23,54,55] , but a large part of this apparent high NUE is related to the externalization of feed production [22,23] . Assigning a NUE coefficient to purchased feed decreased farm NUE and also decreased the differences in farm NUE among intensification levels (Fig. 4). While there is a detailed and well-elaborated methodology for accounting CO2 emissions from the production, transport and use of animal feedstuff [56] , there is no common protocol for assigning a NUE coefficient to off-farm produced feedstuff. We used two coefficients, 50% and 75%, to analyze the sensitivity of farm NUE to accounting for the externalization of feed production, and to indicate the uncertainty in the NUE (and PUE) of purchased feed. Evidently, farm NUE is very sensitive to externalization of feed production; assigning a NUE of 50% to purchased feed nullified the differences among intensification levels in farm NUE. A NUE of 50% is close to the world average NUE in cereal production [57] , but it should be noted that concentrates often contain a large fraction of residues from the food processing industry and only a small percentage of cereals, which makes a comparison with cereal production precarious.
The PUE at farm level was also positively related to the intensity of milk production; mean PUE was 122% for group 1 and 140% for group 10 (Table 1), again with large within-group variability ( Fig. 3(b)). Such high PUE levels and an increasing PUE with intensification level are unlike the situation in most other intensive dairy farms [58,59] . These high PUE values reflect a negative P surplus at farm level and soil P mining, due to P fertilization limits as function of soil P status. Again, the externalization of the cost and (in)efficiency of feed production had a large impact on farm PUE (Fig. 4). Assigning a PUE coefficient to purchased feed strongly decreased overall mean farm PUE and strongly decreased the differences among intensification groups in farm PUE. A relative constancy of P surpluses with intensification level (Fig.  3(a)) has also been reported for dairy farms in Virginia, USA; P input via manure and fertilizers decreased as soil P status increased [60] . However, in many intensive grassland-based dairy systems, both N and P surpluses increase with the intensity of dairy production [23,59,61] .
The decrease in mean P surplus from 1 kg·ha −1 ·yr −1 for group 1 to a mean P deficit of 7 kg·ha −1 ·yr −1 for group 10 reflects the increase in dry matter yield with intensification level (Table S3) and the fact that most intensive farms also had high soil P test values and hence relatively low P application limits. These latter farms are required to lower the soil P status, so as to decrease the risk of soil P leaching and eutrophication of surface waters [62,63] . In the course of time, soil P status will likely decrease and the permissible P input will increase again, without much risk of crop yield losses [64] .
Total NH3 emissions per kg of milk produced decreased with intensification level, from 4.9 for group 1 to 2.7 kg per 1000 kg of milk produced for group 10 (Table 1), a decrease of 0.1 kg per 1000 kg of milk (Table S3). Total NH3 emissions include the emissions from manure stores and following application of manure to farm land. Emissions from applied synthetic N fertilizer (mainly calcium ammonium nitrate) are small ( Table 1). The decrease in total NH3 emissions with intensification level reflects the externalization of the (in)efficiency of animal manure use; NH3 emissions associated with exported manure are not accounted for because these emissions occur on other farms. Quemada et al. [22] indicated that the externalization effects of purchased feed are larger than the externalization effects of manure export. Evidently, further studies are needed to properly account for the externalization of NH3 emissions and the nutrient use effectiveness of exported manure.
The slight increase in N surplus with intensification level most likely reflects the increase in NH3 emissions and other gaseous N losses from animal houses and manure storage systems when total milk production and hence N excretion increases per unit of land area. The N surplus increased on average by 1.5 kg per 1000 kg of milk (Table S3), which translates to 13.2 kg N per dairy cow producing 8800 kg·yr −1 and to about 10% of the annual N excretion of this presumed average dairy cow. This N surplus is considered to be an effectively unavoidable N loss from manure storage, which farmers do not have to account for. These 'unavoidable' N losses are mainly through NH3 emissions, and are less for lowemission houses [65,66] .
Grazing intensity was negatively related to the intensity of milk production ( Table 1). The decline in grazing in intensive grassland-based dairy systems across EU-28 is seen as an unwanted trend [67,68] . In response, FrieslandCampina is now rewarding the milk from dairy farms where the cows graze for ≥ 720 h·yr −1 . Intensification level was also associated with a decrease in the number of young stock per dairy cow, to minimize the feed cost for replacement cattle. It requires farmers to increase the health and longevity of the herd, and also to outsource the rearing of young stock to other farms and countries. We did not account for this outsourcing because of lack of information. Intensification level was also associated with an increase in herbage yield (Table. S3), likely because of reduced grazing losses, and the pressure to maximize farm-grown feed production [69] .
Mean NUE and PUE of the herd significantly increased with intensification level (Table 1, Table 2, and  Table S3), which reflects in part the increase in milk production per cow and the decrease in the number of young stock per cow. The increases in NUE and PUE also reflect the pressure that intensive farms feel to lower the N and P input from purchased feed to the farm, so that the N and P excretion and the need to export manure N and P decrease. Feed suppliers facilitate the optimization of N and P contents of animal feed in response to farm-specific requirements. The mean NUE and PUE of the soil/crop compartment also significantly increased with intensification level, which reflects the effects of N and P application limits and the increase in herbage yield with intensification level. A mean NUE of 66%-71% for grassland and maize fields combined is relatively high compared to the world average NUE of crop production systems [57] .
Other countries in Europe, North America and Oceania with intensive grassland-based dairy production face rather similar challenges with reducing N and P surpluses and with increasing the NUE and PUE of dairy production. Notably Denmark has made progress and decreased N surpluses and increased NUE greatly from mid-1980s in response to the implementation of series of agri-environmental policies [70] . Dairy farms in New York state (USA) were able to decrease N and P surpluses by 30%-50% (per kg of milk produced) between 2004 and 2013 through a combination of incentives and improved nutrient management practices [71][72][73] . Key to such improvements is accurate monitoring of N and P inputs and outputs, and nutrient management advice to farmers.

The KringloopWijzer model as nutrient management instrument
A wide range of nutrient management tools have been developed in dairy production regions of the industrialized world during the last few decades in response to regulatory requirements, market access requirements and the understanding that increasing nutrient use efficiency may go hand in hand with increasing farm profit [74][75][76] . Most of these tools are used as decision support instruments on a voluntary basis. The KringloopWijzer was developed to assess the agronomic and environmental performances of commercial dairy farms in The Netherlands, initially also on a voluntary basis, but since 2016 application of the KringloopWijzer has been mandatory for almost all dairy farms in the Netherlands. This obligation was an initiative of the dairy sector to guide farmers to improve farm management and thereby farm income, as well as to collect information about the performances of individual dairy farms. The obligation has met some objections, especially from dairy farms with a relatively low milk production (< 10,000 kg·ha −1 ·yr −1 ). The data and results of the KringloopWijzer are owned by the farmer and the dairy sector, and are not publicly available.
The basic principles of the various nutrient management tools for dairy farms are rather similar, but there are often differences in accounting for specific inputs and outputs, in assessing farm internal nutrient use efficiencies (feed, herd, manure, soil/crop), and in estimating nutrient losses from the various compartments of dairy farms [74,75] . The KringloopWijzer model is relatively advanced in estimating farm internal nutrient flows and use efficiencies. Most of the input and output data have been standardized, and information from the dairy industry, feed suppliers, fertilizer suppliers, and accounting and registration offices are automatically uploaded. Manure export has to be quantified (weighed), the composition analyzed by a certified laboratory, and these data are automatically registered at the registration office and in the KringloopWijzer. Most uncertain items include stock changes (feed, manure) between years, and exchange of feed stuff among farms and local producers (because these flows are often not certified). Although assessments are made on an annual basis, effects of uncertainties in stock changes fade away over time, because increases in stocks in one year are recorded as inputs in the next year. Partitioning of N and P losses and greenhouse gas emissions are estimated farm-and soil specific, but a in rather simple manner. The KringloopWijzer model is used for decisions at strategic and tactical management levels, often in combination with simulation models like DairyWise [77] and economic/financial management planning.

CONCLUSIONS
A unique database with empirical data from > 5000 dairy farms over three years allowed us to examine the effects of the intensification level of grassland-based dairy production on a range of farm performance indicators. We used data of 2856 farms from the last year of the test phase (2015) because of the large number of farms in this year, and the low number of rejected farm records during the screening tests (because of the learning effect). All farms are required to adopt agri-environmental regulations, which are related in part to the level of intensification of dairy production.
Farms were ordered in ascending milk production per ha of farmland and then divided in 10 equal groups (with 286 farms each). Milk production ranged from a mean of 11 Mg·ha −1 ·yr −1 for group 1 to 30 Mg·ha −1 ·yr −1 for group 10. This wide range allowed us to analyze changes in farm characteristics and performances with the intensification level of dairy production. Additional correlation and regression analyses indicated that the intensity of milk production was related to changes in performance indicators. With an increase in intensification level, feed purchase increased, indicating that dairy production became less land-bound. Due to the regulated manure N and P applications, manure export increased with intensification level. NUE and PUE at farm level increased with intensification level; these increases were to a significant fraction related to the externalization of feed production. Interestingly, N and P surpluses did not change much with intensification level, because of the legal N and P application limits and the obligatory export of manure surpluses.
For long, intensification of dairy production has been associated with increasing nutrient surpluses, decreasing nutrient use efficiencies and accumulation of nutrients in the soil. Our study shows that the intensification level of dairy production is not necessarily related to increasing nutrient surpluses, decreasing nutrient use efficiencies, and increasing soil nutrients stocks, providing the enforcement of N and P application limits, obligatory export of manure N and P surpluses to other farms (and countries), and optimization of whole-farm nutrient management. The first two main requirements have been enforced by the Netherlands government and the last one by the dairy sector itself, in part through the KringloopWijzer model.
The externalization of the (in)efficiencies associated with purchased feed, export of manure, and to a lesser extent the raising of young stock contributed to the apparent increase in NUE and PUE at farm level with increasing milk production intensity. We partly corrected for these externalization effects, but note that there is currently no common protocol for accounting externalization effects. We recommend that such common protocol is made, because intensification of dairy production is a worldwide phenomenon, and farm NUE, PUE, N surplus and P surplus are increasingly seen as important performance indicators.
Intensive livestock production is debated because of its roles in resource use, N and P pollution of the environment, greenhouse gas emissions, and their impacts on landscape, birdlife and biodiversity. Our study contributes to this debate by presenting data and insights from intensive grassland-based dairy farms in The Netherlands that greatly differ in intensity and that are challenged by the government, dairy industry and society at large to improve production performance to be able to obtain and prolong a license to produce.