Article

# Modeling Study on Optimizing Water and Nitrogen Management for Barley in Marginal Soils

Renaldas Žydelis <fnref n="1" />,<fnref n="2" />,* ID, Rafaella Chiarella <fnref n="3" />, Lutz Weihermüller <fnref n="3" /> ID, Michael Herbst <fnref n="3" />, Evelin Loit-Harro <fnref n="4" />, Wieslaw Szulc <fnref n="5" />, Peter Schröder <fnref n="6" /> ID, Virmantas Povilaitis <fnref n="1" />, Michel Mench <fnref n="7" /> ID, Francois Rineau <fnref n="8" /> ID, Eugenija Bakšienė <fnref n="1" /> ID, Jonas Volungevičius <fnref n="1" /> ID, Beata Rutkowska <fnref n="5" /> ID and Arvydas Povilaitis <fnref n="2" /> ID

<figure><img src="image_1.png" /></figure>

### check for updates

Academic Editor: Fucang Zhang

Received: 27 January 2025

Revised: 12 February 2025

Accepted: 20 February 2025

Published: 25 February 2025

**Citation:** Žydelis, R.; Chiarella, R.; Weihermüller, L.; Herbst, M.; Loit-Harro, E.; Szulc, W.; Schröder, P.; Povilaitis, V.; Mench, M.; Rineau, F.; et al. Modeling Study on Optimizing Water and Nitrogen Management for Barley in Marginal Soils. *Plants* **2025**, *14*, 704. https://doi.org/10.3390/plants14050704

**Copyright:** © 2025 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/).

**Abstract:** Water and N availability are key factors limiting crop yield, particularly in marginal soils. This study evaluated the effects of water and N stress on barley grown in marginal soils using field trials and the AgroC model. Experiments from 2020 to 2022 in Lithuania with spring barley cv. KWS Fantex under two N fertilization treatments on sandy soil provided data for model parameterization. The AgroC model simulated barley growth to assess yield potential and yield gaps due to water and N stress. Potential grain yields (assuming no water or N stress) ranged from 4.8 to 6.02 t DW ha⁻<fnref n="1" />, with yield losses up to 54.4% assuming only N stress and 59.2% assuming only water stress, even with the N100 treatment (100 kg N ha⁻<fnref n="1" /> yr⁻<fnref n="1" />). A synthetic case study varying N fertilization from 0 to 200 kg N ha⁻<fnref n="1" /> yr⁻<fnref n="1" /> showed that increasing N still enhanced yield, but the optimal rate of 100–120 kg N ha⁻<fnref n="1" /> yr⁻<fnref n="1" /> depended on climatic conditions, leading to uncertainty in fertilization recommendations. This study underscores the importance of integrating advanced modeling techniques with sustainable agricultural practices to boost yield potential and resilience in marginal soils. Incorporating remote sensing data to capture soil and crop variability is recommended for improving simulation accuracy, contributing to sustainable agriculture strategies in the Baltic–Nordic region.

**Keywords:** barley; crop modeling; marginal soil; nitrogen stress; water stress

## 1. Introduction

Changing future climatic conditions is an urgent global challenge impacting agricultural systems and food production worldwide [1]. It is projected that by 2050 the world's population will surpass 9.5 billion, and consequently, food production needs to increase

<footnote n="1">¹ Institute of Agriculture, Lithuanian Research Centre for Agriculture and Forestry, 58344 Kėdainiai, Lithuania; virmantas.povilaitis@lammc.lt (V.P.); eugenija.baksiene@lammc.lt (E.B.); jonas.volungevicius@lammc.lt (J.V.)</footnote>

<footnote n="2">² Department of Water Engineering, Vytautas Magnus University, 44248 Kaunas, Lithuania; arvydas.povilaitis@vdu.lt</footnote>

<footnote n="3">³ Agrosphere Institute (IBG-3), Forschungszentrum Jülich GmbH, 52428 Jülich, Germany; r.chiarella@fz-juelich.de (R.C.); l.weihermueller@fz-juelich.de (L.W.); m.herbst@fz-juelich.de (M.H.)</footnote>

<footnote n="4">⁴ Field Crops and Plant Biology, Estonian University of Life Sciences, 51006 Tartu, Estonia; evelin.loit-harro@emu.ee</footnote>

<footnote n="5">⁵ Institute of Agriculture, Warsaw University of Life Sciences—SGGW, 02-787 Warsaw, Poland; wieslaw_szulc@sggw.edu.pl (W.S.); beata_rutkowska@sggw.edu.pl (B.R.)</footnote>

<footnote n="6">⁶ Department Experimental Environmental Simulation, Helmholtz Center for Environmental Health, Neuherberg, 85764 Oberschleißheim, Germany</footnote>

<footnote n="7">⁷ University of Bordeaux, INRAE, Biogeco, Bat B2, Allée G. St-Hilaire, F-33615 Pessac CEDEX, France; michel.mench@inrae.fr</footnote>

<footnote n="8">⁸ Environmental Biology, Centre for Environmental Sciences, Hasselt University, 3590 Diepenbeek, Belgium; francois.rineau@uhasselt.be</footnote>

<footnote>* Correspondence: renaldas.zydelis@lammc.lt</footnote>
by 70–85% [2]. Given the scarcity of available arable land within Europe, the effective utilization of previously underutilized or marginal land (except those for biodiversity conservation) becomes a crucial strategy to keep production sufficiently high and even to increase production [3].

Fluctuations in climatic conditions during the growing season directly affect crop physiological processes, including growth, development, and finally, overall productivity in terms of total yield but also yield quality. Crop development mainly depends on air temperature and available water in the root zone, and therefore, fluctuations in these conditions directly affect the length of the growing season and overall growth [4]. Radiation, particularly sunlight, on the other hand, affects the assimilation rate by providing energy for photosynthesis [5], while the amount of precipitation and the water holding capacity of soil determine soil water availability for crop growth and transpiration [6]. Vogel et al. [7] highlighted that climatic factors occurring during the growing season at the global scale can elucidate 20–43% of maize, soybean, rice, and spring wheat yield variation, which is attributed to the explained climate-dependent variance. However, the impact severity of climate change on crop production varies between regions and crops.

Nitrogen (N) is essential for plant growth and development, significantly influencing yield formation patterns [8]. At the same time, N is the most limiting macronutrient for crop production in most agricultural areas worldwide. Consequently, farmers often tend to overfertilize with N to ensure maximum yields and avoid yield reductions associated with N deficiency [9]. In addition, the overall efficient use rate of N fertilizers is usually around 50%, meaning the remaining unused N potentially leaches after crop harvest [10]. Continuous overfertilization of crops without considering plant needs has led to nutrient pollution becoming a major environmental problem for groundwater and surface waterbodies [11].

Barley (*Hordeum vulgare* L.) holds a significant position among cereals, serving an essential role in both animal feeding and human consumption (e.g., in malting) [12]. Barley stands as one of the most ancient, cultivated crops worldwide, and due to its adaptability to various climatic conditions, it can be grown in highly productive high-input agricultural systems but also in subsistence or low-input agriculture, making it a versatile and resilient crop [13]. In Europe, over the past 5 years, the total barley cultivation area remained relatively constant, varying in the range of 22.5–24.2 million ha⁻¹. Hereby, slightly more than half (54.8–57.6%) of the barley cultivation area is concentrated in Eastern Europe, while Northern, Southern, and Western Europe have similar barley cultivation shares, accounting for 13.8–14.4%, 13.8–14.3%, and 15.3–16.3%, respectively [14]. Despite the fact that Eastern Europe accounts for the major part of the barley area, the grain yields in this region are the lowest with 2.6–3.1 t ha⁻¹. A similar yield, i.e., 2.5–4.1 t ha⁻¹, is obtained in Southern Europe, while in Northern and Western Europe, barley yield is higher with approximately 4.3–5.9 and 5.8–6.9 t ha⁻¹ [14]. Northern Europe is generally well suited for barley cultivation due to its favorable climate and fertile soils. However, at present, the main limitations of barley growing in Northern Europe are the occurrence of drought, rain during sowing and harvest, and heat stress [15]. According to the global yield gap atlas (www.yieldgap.org), the barley yield potential in Northern Europe varies between 6.4 and 7.2 t ha⁻¹, while the relative yield gap (difference between actual and yield potential) is 35.3–62.9%, respectively.

Despite barley being a rapidly maturing and relatively drought-tolerant crop, the relative yield gap is likely to be much higher when cultivated on marginal soils compared to cultivation on fertile soils. In this context, marginal soils can be defined by their biological, physical, economic, and environmental limitations or constraints for agricultural production. For instance, it refers to land that does not permit sustainable food production or normal utilization [16]. According to the Muencheberg Soil Quality Rating (SQR) system,
46% of European land is classified as “marginal” [17]. In addition, according to the EUSO soil health dashboard, 61% of EU soils are in an unhealthy state. Based on both sources, the field experiment at the barley growth site described below falls within the poorest area of Lithuanian soils, where sandy soil textures prevail. As a consequence, these soils have low water and nutrient retention and are also characterized as acidic soils. Therefore, as the global population grows and pressure on productive agricultural land increases, there is a demand to bring back more non-fertile soils into agricultural use. Despite the favorable conditions for barley cultivation under Baltic-Nordic environments, there is still a lack of experimental data, especially regarding yield potential and yield gaps when barley is grown under marginal soil conditions. The results of such a study could be highly beneficial in providing guidance for barley breeding tailored to low-fertility soil conditions.

Obviously crop growth, biomass, and yield development are significantly influenced by prevailing weather and soil conditions, as well as by management practices, such as sowing date, seeding density, and fertilizer application [18]. Under field conditions, identifying these limiting factors (such as the interaction between N and water stress) on plant growth and yield is challenging, especially when analyzing barley production in marginal soils. However, the complex relationships among these factors are well captured in state-of-the-art process-based models, and therefore, such models can serve as vital decision-support tools for farmers and farm advisors, enabling the calculation of potential crop yields and yield gaps under current or even future climate conditions [19,20] but also under various management strategies, such as different fertilization levels.

As outlined, there is still a lack of studies analyzing the yield potential and yield gaps of barley production on marginal soils under Baltic-Nordic conditions. To address these issues, the process-based agroecosystem model AgroC [21,22] based on a multiyear barley experiment with contrasting climatic conditions was calibrated and validated, and subsequently used to simulate barley growth responses to various N fertilization scenarios. The AgroC model was selected for its physically based soil water module, which may offer superior performance compared to traditional water bucket models commonly used in crop growth simulations, particularly for modeling barley growth under nitrogen and water stress conditions

The objectives of this study were as follows: (1) calibrate and evaluate the performance of the AgroC model using three-years' field experimental data from barley cultivation on a marginal soil in the Nemoral climate zone; (2) to assess the potential spring barley yield and to untangle and quantify the confounding impacts of N and water stress, and (3) to unravel the impacts of increasing N application rates on barley yield, N uptake, nitrate leaching, and the intensity of abiotic stress.

## 2. Results and Discussion

### 2.1. Field Observations

The barley growing period lasted 108 days in 2020, 96 days in 2021, and 98 days in 2022. The total aboveground barley biomass (TAB) yield during this period varied significantly, from 4.15 ± 0.49 t ha⁻¹ in 2020 to 4.97 ± 0.68 t ha⁻¹ in 2021 for N100 treatment plots. In contrast, the TAB yield for the N0 (non-fertilized) treatment plots was approximately half of that found for the N100 plots, ranging from 1.82 ± 0.21 t ha⁻¹ in 2020 to 2.49 ± 0.70 t ha⁻¹ in 2021. These results indicate substantial differences between the growing seasons and fertilization conditions. From sowing to harvest, the mean air temperatures were 15.9 °C in 2020, 18.6 °C in 2021, and 17.9 °C in 2022 (hence 0.3 °C lower, 1.9 °C higher, and 1.7 °C higher compared to the 1990–2020 long-term average, respectively). Precipitation amounts during the barley growing seasons were also very variable. In 2020, the total precipitation amount was only 225.4 mm (94.7% of the long-term average), while in 2021 and 2022,
similar amounts were recorded, with 327.8 (147.7% of the long-term average) and 346.2 mm (145.5% of the long-term average), respectively. Consequently, the 2020 season will be referred to as "normal" for both temperature and precipitation regimes, while 2021 and 2022 are referred to as "hot" and "wet".

According to the official statistics in Lithuania [23], the average yield of spring barley across the country was 4.22, 3.30, and 3.81 t ha⁻¹ in 2020, 2021, and 2022, respectively. These data were used to estimate the difference between the barley yields obtained on farmers' fields and those from the field experiments with the N100 fertilization rate. The barley yields from the field experiments were lower than the average yield achieved across the country by 50, 33, and 41% in 2020, 2021, and 2022, respectively. These differences are mainly caused by the marginal soils the barley grew on in the experiments. In contrast, when comparing the yields from the field experiments with those from the same area (e.g., around Vilnius region), where predominantly similar soil types (Arenosols) were reported, the yield was only 20–30% lower in 2020–2022. This difference was likely the result of slightly lower fertilization rates used in our experiment compared to those typically used by farmers in the region.

## 2.2. Model Calibration and Validation

### Soil Water Content

The measured soil water content (SWC) across five depths during the three year experimental period and for the N100 plot is illustrated in Figure 1. Additionally, statistical indicators, such as d, RMSE, and BIAS, are detailed in Table 1. Changes in measured SWC were higher in the arable Ap (0–30 cm) and B1 (30–50 cm) soil horizon compared to deeper horizons, such as B2 (50–78 cm), 2C(k) (78–105 cm), and 2C (105–120 cm). During the barley growing season, SWC values at 15 cm depth varied within a range of 0.05–0.29 cm³ cm⁻³, while at 40 cm depth, SWCs varied from 0.06 to 0.18 cm³ cm⁻³. At a larger depth, the SWC ranged even more narrowly, and the soil profile was never fully saturated or even near saturation. However, in the year 2020, which was referred to as "normal", the SWC values across all soil horizons were lower and exhibited less fluctuation in relation to the hot/wet years of 2021–2022. Despite the varying trends in SWC fluctuations, the AgroC model matched the SWC dynamics well after calibration for the experimental years. Slightly better statistical outcomes were achieved during the calibration period (2020 and 2022) compared to the validation period (2021). The agreement between simulated and measured SWC at a 15 cm depth was reasonably good and similar for both the calibration and validation periods, with d values of 0.620 and 0.788, RMSE of 0.051 and 0.067 cm³ cm⁻³, and BIAS of 0.066 and 0.040, respectively. At the B1 horizon at 40 cm depth, the correspondence between measured and simulated SWC was acceptable, showing slightly better performance during the calibration period (d = 0.542, RMSE = 0.031 cm³ cm⁻³, BIAS = -0.001) compared to the validation period (d = 0.391, RMSE = 0.039 cm³ cm⁻³, BIAS = -0.006). In contrast, at 90 cm depth, the SWC values were underestimated for both periods (calibration and validation), particularly for 2021 and 2022 (Figure 1). At the deepest measured depth for SWC (110 cm), which corresponded to the lowest rooting depth of the barley, the SWC values varied the least, and the agreement between measured and simulated values was fairly good as well. For the calibration period, the average d index was 0.426, with a small RMSE of 0.015 cm³ cm⁻³ and a BIAS close to zero at 0.011. For the validation period, the corresponding values were similar with a d index of 0.455, RMSE of 0.008 cm³ cm⁻³, and a BIAS of -0.003.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Comparison between observed (dots) and simulated (lines) volumetric soil water content, θ (cm³ cm⁻³) at 15, 40, 60, 90, and 110 cm depths for the N100 treatment during calibration (2020 and 2022) and validation (2021) periods for rainfed barley. SAT is the saturated water content, PWP is the water content at the permanent wilting point, and SAT-PWP is the plant available water (PAW) (all in cm³ cm⁻³). Error bars indicate standard deviation of the measurements.</figcaption></figure>

<table><caption>**Table 1.** Statistical values between the measured vs. simulated data for the calibration (2020 and 2022) and validation period for rainfed barley (treatments: N0 and N100) at the experimental Voké site.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">d</th><th colspan="2">RMSE</th><th colspan="2">BIAS</th></tr><tr><th>N0</th><th>N100</th><th>N0</th><th>N100</th><th>N0</th><th>N100</th></tr></thead><tbody><tr><td colspan="7">2020</td></tr><tr><td>Leaf area index</td><td>0.704</td><td>0.834</td><td>0.663</td><td>0.535</td><td>-0.619</td><td>0.286</td></tr><tr><td>TAB (t ha<sup>-1</sup>)</td><td>0.982</td><td>0.980</td><td>0.128</td><td>0.422</td><td>-0.029</td><td>-0.321</td></tr><tr><td>Leaf (t ha<sup>-1</sup>)</td><td>0.702</td><td>0.513</td><td>0.193</td><td>0.444</td><td>0.074</td><td>0.347</td></tr><tr><td>Stem (t ha<sup>-1</sup>)</td><td>0.892</td><td>0.587</td><td>0.172</td><td>0.446</td><td>-0.163</td><td>-0.420</td></tr><tr><td>Storage organs and GY (t ha<sup>-1</sup>)</td><td>0.725</td><td>0.928</td><td>0.241</td><td>0.315</td><td>0.258</td><td>0.242</td></tr><tr><td>SWC at 15 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>—</td><td>0.788</td><td>—</td><td>0.030</td><td>—</td><td>0.040</td></tr><tr><td>SWC at 40 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>—</td><td>0.775</td><td>—</td><td>0.019</td><td>—</td><td>-0.017</td></tr><tr><td>SWC at 60 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>—</td><td>0.907</td><td>—</td><td>0.004</td><td>—</td><td>0.002</td></tr><tr><td>SWC at 90 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>—</td><td>0.159</td><td>—</td><td>0.013</td><td>—</td><td>0.012</td></tr><tr><td>SWC at 110 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>—</td><td>0.101</td><td>—</td><td>0.019</td><td>—</td><td>0.018</td></tr></tbody></table>
<table><caption>**Table 1.** *Cont.*</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">d</th><th colspan="2">RMSE</th><th colspan="2">BIAS</th></tr><tr><th>N0</th><th>N100</th><th>N0</th><th>N100</th><th>N0</th><th>N100</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">2021</td></tr><tr><td>Leaf area index</td><td>0.719</td><td>0.806</td><td>0.795</td><td>0.875</td><td>-0.787</td><td>0.037</td></tr><tr><td>TAB (t ha<sup>-1</sup>)</td><td>0.998</td><td>0.925</td><td>0.051</td><td>1.197</td><td>0.034</td><td>-0.973</td></tr><tr><td>Leaf (t ha<sup>-1</sup>)</td><td>0.728</td><td>0.755</td><td>0.196</td><td>0.365</td><td>0.169</td><td>0.323</td></tr><tr><td>Stem (t ha<sup>-1</sup>)</td><td>0.802</td><td>0.570</td><td>0.159</td><td>0.938</td><td>-0.100</td><td>-0.786</td></tr><tr><td>Storage organs and GY (t ha<sup>-1</sup>)</td><td>0.975</td><td>0.643</td><td>0.166</td><td>0.624</td><td>0.027</td><td>-0.430</td></tr><tr><td>SWC at 15 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.788</td><td>-</td><td>0.067</td><td>-</td><td>0.040</td></tr><tr><td>SWC at 40 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.391</td><td>-</td><td>0.039</td><td>-</td><td>-0.006</td></tr><tr><td>SWC at 60 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.468</td><td>-</td><td>0.019</td><td>-</td><td>-0.005</td></tr><tr><td>SWC at 90 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.339</td><td>-</td><td>0.031</td><td>-</td><td>-0.025</td></tr><tr><td>SWC at 110 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.455</td><td>-</td><td>0.008</td><td>-</td><td>-0.003</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">2022</td></tr><tr><td>Leaf area index</td><td>0.882</td><td>0.949</td><td>0.384</td><td>0.429</td><td>-0.349</td><td>0.295</td></tr><tr><td>TAB (t ha<sup>-1</sup>)</td><td>0.629</td><td>0.972</td><td>0.772</td><td>0.527</td><td>0.656</td><td>0.485</td></tr><tr><td>Leaf (t ha<sup>-1</sup>)</td><td>0.389</td><td>0.582</td><td>0.334</td><td>0.423</td><td>0.314</td><td>0.385</td></tr><tr><td>Stem (t ha<sup>-1</sup>)</td><td>0.438</td><td>0.992</td><td>0.248</td><td>0.113</td><td>0.248</td><td>0.069</td></tr><tr><td>Storage organs and GY (t ha<sup>-1</sup>)</td><td>0.408</td><td>0.914</td><td>0.411</td><td>0.470</td><td>0.266</td><td>0.057</td></tr><tr><td>SWC at 15 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.452</td><td>-</td><td>0.071</td><td>-</td><td>0.093</td></tr><tr><td>SWC at 40 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.308</td><td>-</td><td>0.043</td><td>-</td><td>-0.006</td></tr><tr><td>SWC at 60 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.293</td><td>-</td><td>0.035</td><td>-</td><td>-0.027</td></tr><tr><td>SWC at 90 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.108</td><td>-</td><td>0.072</td><td>-</td><td>-0.065</td></tr><tr><td>SWC at 110 cm (cm<sup>3</sup> cm<sup>-3</sup>)</td><td>-</td><td>0.751</td><td>-</td><td>0.011</td><td>-</td><td>0.004</td></tr></tbody></table>

Note: SWC—soil volumetric water content, *d*—index of agreement, RMSE—root mean squared error, BIAS—systematic bias index.

## 2.3. Leaf Area Index and Barley Development Stages

The comparisons between simulated and measured barley development stages (DVS) after calibrating the phenology are presented in Figure 2. Throughout the years 2020–2022, crop germination, aided by favorable weather conditions, typically spanned 7–10 days until emergence. The barley flowering (BBCH = 60–69 or DVS = 1) and physiological maturity (BBCH = 92–99 or DVS = 2) stages were recorded 60–68 and 95–107 days after seeding, respectively. The corresponding simulated values for the 2020–2022 period for flowering and physiological maturity varied between 56 and 68 and between 92 and 104 days after seeding, indicating a strong concordance between the simulated and measured DVS. During the barley growing seasons from 2020 to 2022, only a few days of crop phenology differences were observed between the N100 and N0 treatments. Crop phenology is primarily determined by genetic characteristics and environmental factors, such as temperature, rainfall, and photoperiod [24]. However, the slight differences in crop development observed over those few days were likely influenced by the varying nutrient conditions.

During the model calibration (2020 and 2022) and validation (2021) periods, there was a close match between the simulated and observed crop leaf area index (LAI) (Figure 2, Table 1). During calibration, the simulated LAI values for the N100 treatment were slightly overestimated, as indicated by a *d* index of 0.892, an RMSE of 0.482 m² m⁻², and a BIAS of 0.291. However, for the validation period, the simulated results were closely aligned with the measured values, demonstrating a high *d* index of 0.925, a RMSE of 0.875 m² m⁻², and a BIAS close to 0 at 0.037. For the N0 treatment, the simulated LAI values were slightly underestimated but fell within the standard deviation limits. During the calibration period, the *d* index was 0.793, with a RMSE of 0.524 m² m⁻² and a BIAS of -0.484. Similar LAI statistical results were observed for the validation period, with a *d* index of 0.719, a RMSE
of 0.795 m² m⁻², and a BIAS of -0.787. Here, it has to be noted that the mismatch in model prediction for the individual years were likely caused by different plot locations within the field as the soils in the field were not fully homogeneous and even small differences in soil texture or bulk density will have a large impact in the soil hydraulic characteristics of the extreme coarse and well sorted sandy soil.

<figure><img src="image_3.png" /><figcaption>**Figure 2. (top)** Observed and simulated dry matter biomass of individual vegetative and reproductive parts of barley (leaves, stems, and grain yield). (middle) Seasonal variation in the green leaf area index (LAI) and (bottom) barley development stages (DVS) across 3 growing seasons, 2020–2022. In all figures, observed values are represented by dots, with error bars indicating the standard deviation, while simulated values are shown as lines for the N100 (left panel) and N0 (right panel) fertilization.</figcaption></figure>

## 2.4. Partitioning of Total Above-Ground Biomass, Leaf, Stem, and Grain Yield

The measured and simulated total aboveground biomass (TAB) as well as biomass of the individual organs for two contrasting N treatments, N0 and N100, are presented in Figure 2 and Table 1. The statistical measures showed that for the N100 treatment, the model slightly better reproduced seasonal variation of the TAB and its individual components during the calibration period compared to the validation period. For N0 treatment, the results were significantly better during the validation period, with d values very close to 1.

For example, for the N100 treatment, during the calibration period, the d value was high at 0.976, associated with a low RMSE of 0.475 t ha⁻¹ and BIAS = 0.082 for TAB, but the d value (=0.925) decreased for the validation period in terms of RMSE = 1.197 t ha⁻¹ and BIAS = -0.973. Analyzing the grain yields for the N100 treatment, the trends were similar to the observations for TAB. However, during the validation period, there was a slight overestimation (see Table 1) of biomass and grain yield (GY). These differences in simulated TAB and GY for the validation period likely resulted from changes in grain components,
specifically the grain number and 1000 grain weight. In 2021, there was a larger grain count but lower grain weight.

Specifically, during the calibration period, the weight of 1000 grains ranged from 40 to 47 g, whereas during the validation period, it was only 30–35 g. Additionally, variations in the grain nutrient content were observed. During the calibration period, the N content in barley grain was approximately 2.36%, compared to a lower N content of about 2.08% during the validation period. Despite the high N100 fertilization rate, this indicated greater N stress in 2021 (validation period) compared to the calibration periods of 2020 and 2022. On the other hand, the P and K contents in both stems and leaves remained similar throughout the 2020–2022 period.

As indicated by the statistical indices listed in Table 1, the AgroC model showed varying degrees of accuracy in simulating leaf and stem biomass for the N100 treatment. For example, for stem biomass, better statistical indices were obtained during the calibration period with a d index of 0.790, RMSE of $0.570 \ t \ ha^{-1}$, and BIAS of -0.176, compared to the validation period with a d index of 0.570, RMSE of $0.928 \ t \ ha^{-1}$, and BIAS of -0.786. Conversely, for leaf biomass simulation, the model showed better predictions during the validation period compared to the calibration period.

The AgroC model performed exceptionally well in simulating TAB for the N0 treatment during the validation period, with very good statistical measures. This contrasted with the calibration period, in which the performance was not as strong, particularly in 2022. During the calibration period, the d index was 0.806, RMSE was $0.450 \ t \ ha^{-1}$, and BIAS was 0.314. In comparison, the validation period showed a high d index of 0.998, a very low RMSE of $0.051 \ t \ ha^{-1}$, and BIAS of 0.034, indicating almost perfect model accuracy. Similar trends were observed when analyzing the dynamics of storage organs for the N0 treatment. The model performed better during the validation period compared to the calibration period. During the calibration period, the model simulated storage organs with moderate agreement, a d index of 0.567, RMSE of $0.326 \ t \ ha^{-1}$, and BIAS of 0.262. In the validation period, these measures improved significantly to a d index of 0.975, RMSE of $0.166 \ t \ ha^{-1}$, and BIAS of 0.027. Conversely, the AgroC model performed better in simulating leaf and stem biomass in unfertilized plots during the calibration period compared to the validation period. There was a noticeable increase in RMSE and a shift from overestimation to underestimation during the validation period.

Overall, the accuracy of the AgroC model during calibration and validation should be discussed with the understanding that this model was designed to predict crop yields at the field and field plot scale (point simulation). As such, the field was assumed to be uniform, without accounting for spatial differences in crop development caused by small scale changes in soil characteristics and/or crop management. Calibration and validation results for both N0 and N100 treatments indicated that while the model was highly effective, some underestimations and overestimations of biomass and individual organs may be due to changes in N availability under different environmental conditions. Moreover, sandy soils (where the field experiments were conducted) may exhibit significant morphological, chemical, and physical variability at various scales [25] even though they look fairly homogeneous on the other side. Previous studies have shown that soil parameters, such as soil organic carbon (SOC), but also soil texture, such as clay and sand content, can vary by more than twofold even within a small 3–4 ha field [26]. Therefore, the spatial variability of soil parameters could likely influence the model results in our study as the plots used in the three years were not the same. These factors could contribute to discrepancies between measured and simulated values, emphasizing the importance of considering soil heterogeneity and environmental variability in future model predictions even for apparently homogeneous sandy soils used in this study.
## 2.5. Estimated Water Stress

Despite barley being considered relatively tolerant to abiotic stress compared to other cereal crops [27], specifically water stress might be an important growth limitation as the field experiments in this study were performed on a highly sandy Arenosol soil (refer to Supplementary Figure S1) with a low water retention capacity. Therefore, the simulated water stress over the three contrasting growing seasons was analyzed and is presented in Figure 3. The average water availability over the entire root zone, $\alpha_{avg}$ (-), which ranges from 1 to 0, was used to identify periods of water stress. $\alpha_{avg}$ is thereby the ratio of actual to potential transpiration and values of 1 indicate no water stress, while smaller values indicate water stress. Based on the plots for the years 2020–2022, the occurrence of water stress and water availability reflected the precipitation and potential evapotranspiration distribution over the growing seasons. The cumulative potential evapotranspiration from 2020 to 2022 was approximately 390–400 mm during the growing season, while the precipitation amount varied between 225.4 and 366.5 mm, indicating an insufficient water supply for barley growth. Despite the different fertilization treatments (N0 and N100) and associated differences in biomass production, the water stress experienced was nearly identical for both treatments. For example, in 2020, the average water stress $\alpha_{avg}$ for the entire barley growing season was 0.76 for the N0 treatment and only slightly smaller with 0.77 for the N100 treatment. In contrast, in 2021, $\alpha_{avg}$ showed higher values for both treatments with 0.86 for N0 and 0.84 for N100, indicating lower water stress compared to 2020. In 2022, $\alpha_{avg}$ values were 0.75 for N0 and 0.74 for N100, indicating the largest stress across all 3 years.

When analyzing specific barley growth periods separately, namely the vegetative stage (from sowing to flowering) and the reproductive stage (from flowering to physiological maturity), it was observed that the average water stress during the vegetative period was fairly consistent across all years. The $\alpha_{avg}$ values during the vegetative period were as follows: in 2020, N0—0.86 and N100—0.85; in 2021, N0—0.80 and N100—0.80; and in 2022, N0—0.86 and N100—0.86. However, the reproductive period experienced more severe water stress in 2020 ($\alpha_{avg}$ values of 0.63 and 0.64 for N0 and N100, respectively) and in 2022 ($\alpha_{avg}$ values of 0.48 and 0.52 for N0 and N100, respectively). In contrast, the reproductive period in 2021 showed much less water stress, with average $\alpha_{avg}$ values ranging from 0.91 to 0.95 for both treatments. It is also important to note that the intensity of $\alpha_{avg}$ reached its peak at different times in different years. For example, in 2020, the highest $\alpha_{avg}$ was simulated during multiple periods, reaching severe stress conditions with an $\alpha_{avg}$ of 0.2–0.3 during the stem elongation–booting growth stage (BBCH31–BBCH47) and even dropping to 0.08 during the flowering–development of grain growth stage (BBCH61–BBCH71). In 2021, the highest $\alpha_{avg}$ was simulated only during the leaf development stage, up to BBCH13, where water stress ranged between 0.36 and 0.44. In 2022, the highest $\alpha_{avg}$ occurred during the stem elongation stage (0.21–0.35) and the flowering–development of grain stages (0.17–0.36).

According to the FAO [14], barley's water demand per growing period typically ranges between 400 and 650 mm, with the exact amount largely depending on the length of the growing season, environmental conditions, soil type, and the genetic characteristics of the variety. During the germination and early growth stages (BBCH stages), barley requires approximately 50–100 mm of water, accounting for about 15.4% of the total seasonal water requirement. As the crop enters the vegetative period, from tillering to stem elongation, water demand increases, requiring an additional 100–200 mm, which represents around 30.8% of the total water need. The most critical phase is the reproductive period, particularly from heading to grain filling, where the water requirement peaks at 200–250 mm, making up about 38.4% of the total. Finally, during the maturation period, barley needs approximately 50–100 mm, accounting for the remaining 15.4% of the total water demand. Throughout the
growing seasons, under rainfed conditions, the crops consistently experienced water stress due to insufficient water supply in our experiments. However, the intensity of this stress varied from mild to severe. As a result, farmers cultivating crops on marginal soils should consider the potential benefits of irrigation to reduce crop stress and increase economic yields. Currently, irrigation is not widely adopted among farmers in the region, with the irrigation area covering only approximately 6600 hectares, relying primarily on surface water (from ponds) and groundwater sources [23].

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Estimated water availability index ($\alpha_{avg}$) for both fertilization treatments (N0 and N100), along with the difference in $\alpha_{avg}$ between N0 and N100, potential evapotranspiration (ETp), and precipitation (P) for the barley growing seasons from 2020 to 2022.</figcaption></figure>

## 2.6. Estimated N Stress, Uptake, and Nitrate Leaching

Determining the optimal N rate for barley to maximize crop productivity while minimizing environmental risk requires an understanding of barley's N use mechanisms. Therefore, two contrasting N treatments were analyzed in this study: one without fertilization (N0) and one with a locally assumed full fertilization of 100 kg ha⁻¹ (N100) applied nitrogen. Consequently, the simulated N availability index, N uptake, and nitrate leaching out of the soil profile over the three contrasting growing seasons were analyzed and are presented in Figure 4. The average N availability index, which ranged from 1 to 0, was used to identify periods of barley N stress. As for the water stress, values of 1 indicate no N stress, while smaller values indicate increasing N stress. Unlike the analysis of water stress, the estimated N stress differed significantly between the N0 and N100 treatments as expected. In 2020, the average N stress for the entire barley growing season was 0.65 for N0, indicating mild N stress and 0.97 for N100, showing almost no N limitation under
N100 fertilization. In 2021 and 2022, N stress in the N0 plots was higher than in the first year with moderate N stress of 0.54 and 0.51, respectively, while N100 also showed some mild N stress ranging from 0.79 to 0.83. When analyzing specific barley growth phases, no N stress was observed during the vegetative period for N100, except in 2022, when it reached 0.92. During the reproductive period, N stress for N100 was present in all years with 0.94 in 2020, 0.47 in 2021, and 0.65 in 2022, indicating that the second year was that with the largest N stress. In contrast, the N0 treatment showed N stress throughout the entire season, being slightly lower during the vegetative period (0.82 in 2020, 0.72 in 2021, and 0.61 in 2022) and severe during the reproductive period (0.42 in 2020, 0.26 in 2021, and 0.32 in 2022). Again, the second year was the growing season with the largest N stress during the reproductive stage.

The ability of the AgroC model to reproduce the barley N uptake by the grain, stem, and leaf is shown in Figure 4. As illustrated, simulated barley N uptake varied across years and fertilization conditions. During the 2020–2022 period, at barley physiological maturity, observed N uptake in the grain ranged from 51.2 to 58.9 kg ha⁻¹ N for the N100 treatment and was significantly less for the N0 treatment, ranging from 12.8 to 18.9 kg ha⁻¹ N (25 and 32% of N100 treatment). At its peak, N uptake by leaves varied between 13.5 and 20.4 kg ha⁻¹ N for the N100 treatment and was approximately half that amount for the N0 treatment, ranging from 5.1 to 7.4 kg ha⁻¹ N (38 and 36% of the N100 treatment). N uptake by stems in 2020–2022 ranged from 9.1 to 14.1 kg ha⁻¹ N for the N100 treatment and from 2.1 to 4.2 kg ha⁻¹ N for the N0 treatment (23 and 30% of the N100 treatment). In nearly all cases, the AgroC model successfully represented the observed N uptake data. However, the model underestimated N uptake by grains for both treatments in the last year. This discrepancy was likely due to the varying timing of water and N stress during the barley growing seasons (2020–2022), as well as instances when both stressors affected the barley simultaneously.

Differences in the simulated mass of nitrate leached below a 120 cm soil depth were primarily driven by year-to-year climatic variability, as shown in Figure 3, initial N status, and differences in fertilizer treatments. Simulated nitrate leaching was consistently higher in the N100 treatment than in the N0 one across all growing seasons and years. The highest nitrate leaching generally occurred during the winter months, when downward water flow occurs due to precipitation exceeding evaporation and a lack of N uptake by plants. For example, in the N100 treatment, the nitrate concentration in leachate during the 2020–2021 and 2021–2022 winter periods reached nearly 30 mg L⁻¹, while in the N0 treatment, the nitrate concentration in leachate peaked at approximately 20 mg L⁻¹ during the 2020–2021 winter but decreased to 15 mg L⁻¹ in the following winter. This higher nitrate concentration in leachate during the first winter season was likely influenced by the relatively high initial soil inorganic N reserves before barley sowing, which were 71.3 kg N ha⁻¹ (NO3-N + NH4-N) even though the model has been spun up over 5 years to avoid such artifacts. Leaching during the autumn period was also substantial. For example, in the N100 treatment, the N concentration in leachate in 2020–2021 ranged from 9.4 to 21.7 mg L⁻¹, while in the following autumn season, it decreased to 9.3–13.4 mg L⁻¹. In comparison, the N0 treatment showed lower values, ranging from 7.0 to 16.2 mg L⁻¹ in 2020–2021 and 6.7 to 9.3 mg L⁻¹ in the next autumn.

During the summer, nitrate leachate concentrations were lower, ranging from 2.4 to 10.0 mg L⁻¹ for both treatments across all years. In spring, the simulated nitrate concentration in leachate was generally low, except in 2021, when it reached 30.0 mg L⁻¹ for the N100 treatment and 20.7 mg L⁻¹ for the N0 treatment. This higher nitrate leaching in spring 2021 was likely due to the unfavorable environmental conditions in 2020, which
resulted in barley not fully utilizing the applied nitrogen. Consequently, in spring, with the soil not covered by any vegetation, much of the remaining nitrate was leached downwards.

<figure><img src="image_5.png" /></figure>

**Figure 4.** Simulated nitrogen availability index for both fertilization treatments (N100 left panel and N0 right panel), along with seasonal N uptake per organ and nitrate leaching during the barley growing seasons of 2020–2022. In the top and middle figures, observed values are represented by dots with error bars indicating the standard deviation, while simulated values are shown as lines.

## 2.7. Model Application to Assess the Environmental Impact and Effect of Increasing Fertilizer on Barley Yields

### 2.7.1. Yield Potential, Water-Limited Yield, Nitrogen-Limited Yield, N Uptake, and N Leaching

The predicted potential barley yield (without water and N stress) was (in t ha⁻¹) 4.80 in 2020, 6.02 in 2021, and 4.94 in 2022 (Figure 5). These potential yield fluctuations were mainly influenced by the different meteorological conditions between the growing seasons and associated changes in solar radiation and temperature.

In the next step, the simulated yield for the N100 treatment over three growing seasons was compared to potential yield without water stress. The grain yield losses due to water shortage were significant in 2020, reaching 2.61 t ha⁻¹ (−54.4%). In 2021, the yield losses due to water stress were only 0.65 t ha⁻¹ (−10.8%), indicating more favorable meteorological conditions compared to the other years. In 2022, the yield losses amounted to 2.23 t ha⁻¹ (−45.1%), which was slightly smaller than those estimated for the first year.

Nitrogen stress reduced the grain yield significantly by 2.60 t ha⁻¹ (−54.3%) in 2020, by 2.92 t ha⁻¹ (−59.2%) in 2021, and by 3.00 t ha⁻¹ (−49.9%) in 2022, whereby year to year fluctuations were less pronounced compared to water stress impacts. The calculated barley
yield gap on the marginal sandy soil, influenced by the combined effects of water and N stress, was 56.9% in 2020, 50.5% in 2021, and 67.4% in 2022. This indicated that even under favorable environmental conditions (second year with water yield losses of only −10.8%), the grain yield reduction due to the combined stresses was large.

The potential yield and yield gap results provided key insights into the interaction between environmental conditions, N management, and water availability, which are all crucial for optimizing barley production on sandy soils. This study evidenced potential barley grain yields ranging from 4.8 to 6.02 t ha⁻¹ (dry weight) on sandy soil at a latitude of 54° N in the Baltic region. According to the Global Yield Gap Atlas (www.yieldgap.org), the estimated barley yield potential in Lithuania ranged from 7.0 to 7.2 t ha⁻¹, depending on climate zones. However, the Atlas information does not account for soil types, making direct comparisons with our results less applicable. Despite barley's inherent tolerance to abiotic stress [28], our findings showed that N and water stress during growth largely reduced yields on sandy soils across all experimental years. Yield losses were even more severe in unfertilized plots, with losses of 81.2, 83.3, and 81.7% in 2020, 2021, and 2022, respectively, under full N (N0) and water stress. According to the Global Yield Gap Atlas, average barley yield losses in Lithuania are around 49.2%. Across Europe, under highly contrasting climate conditions, rainfed barley absolute yield gaps vary between 1.0 and 7.3 t ha⁻¹. This corresponds to a relative yield gap ranging from 12% to 75% of the potential yield [29]. In general, yield gaps can vary between years [30] as also shown in our study and depend on N and water management [31].

In a next step, the optimal N fertilization rates in which N does not limit the barley yield potential were analyzed using the results from the synthetic modeling study, which showed that the N fertilization rates were increased in steps of 10 kg N ha⁻¹ from no fertilization (N0) to 200 kg N ha⁻¹ (N200). Here, the fertilization was always split into two applications as done for the field experiments. Again, a five-year spin-up prior to the year 2020 was added to account for changes in soil N status due to the different fertilization rates. In Figure 5, the grain yield at harvest, N uptake for the different fertilization levels, as well as the N leaching for the N0, N100, and N200 scenario and individual years are depicted. In 2020, the N fertilization rate for optimal harvest was 110 kg ha⁻¹, which was a slightly higher fertilization amount than used in the field experiments, which was assumed to be optimal for this site. However, in 2021 and 2022 the optimal rate increased to 170 kg ha⁻¹, indicating that higher N inputs were needed under more favourable climatic conditions for optimal barley production. Irrespectively of those high optimal N fertilization rates obtained from the simulation results, farmers should consider the economic and environmental optimal N rate, which are often lower than the agronomic maximum [32]. The simulated barley N uptake by grains showed also variation among the years, likely due to differences in simulated grain yields (Figure 5). In 2020, N uptake under the N0 treatment was 15.9 kg N ha⁻¹, increasing to 54.79 kg N ha⁻¹ with higher N fertilization (up to N110), where it reached saturation. Conversely, in 2021, likely due to favourable soil moisture conditions, N uptake increased from 23.9 to 126 kg N ha⁻¹ for the N0 to the N180 rate. In 2022, N uptake ranged from 9 to 65.4 kg N ha⁻¹ for the same rates. Accordingly, the N fertilization limit was N110 in 2020, while it was much higher in 2021 and 2022, reaching N170. This outcome is highly relevant to Lithuania and other Baltic-Nordic countries, where the political focus is currently on reducing N application to protect surface and groundwater bodies as requested by EU policies [33]. Indeed, farmers often use 'excessive' fertilization as a form of insurance to generate optimal yields, but overfertilization can result from overestimating expected yields or neglecting local soil conditions, which are typically the basis for N fertilizer recommendations. For example, a meta-analysis of N fertilizer experiments in Finland showed that current fertilization
recommendations for spring cereals, including barley, which are based on growers' yield expectations, can lead to significant errors in N management practices [34]. According to Lithuania's national fertilization recommendations, N fertilizer applications to achieve a barley grain yield of 1 t ha⁻¹ are 21 kg of N, 10 kg of P, and 21 kg of K [35]. However, our experiments showed that barley's N requirements were slightly higher, at 23 to 24 kg of N per ton of grain yield. Additionally, these national recommendations are based on a standard average barley yield of 4.4 t ha⁻¹, which are higher than those recorded over this 3-year study years on the marginal soil. Looking at the negative site effects of fertilization in terms of N leaching to groundwater, one can see that the barley that received the highest rate of N fertilization (200 kg ha⁻¹) exhibited the highest NO₃⁻-leaching rate over three growing seasons (Figure 5). As shown, during January to March 2021, nitrate leaching fluctuated between 48 and 56 mg N L⁻¹, exceeding the maximum allowable limit of 50 mg N L⁻¹ [36]. In contrast, nitrate leaching for N100 was approximately half of that, around 30 mg N L⁻¹, while for N0, it was even lower at 20 mg N L⁻¹. In February 2022, nitrate leaching for N200 remained elevated, ranging from 34 to 44 mg N L⁻¹, and again exceeded the allowable limit in March, reaching 56 mg N L⁻¹. The risk of nitrate leaching is closely linked to the excessive use of fertilizers, particularly occurring after the growing season and within the crop cycle [37].

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Simulated barley grain yield potential (no stress), water-limited, nitrogen-limited (top row), and simulated nitrogen uptake by grains (middle row) and simulated nitrate leaching among three contrasting treatments: N0, N100, and N200.</figcaption></figure>
### 2.7.2. Estimated Nitrogen Stress and Nitrogen Use Efficiencies

The average N stress intensity over the entire barley growing season along with various indicators reflecting the N use efficiency are plotted for the different N fertilization rates in Figure 6. As shown, N stress intensity was influenced by both fertilization rates and yearly environmental conditions and reflected exactly what has been discussed for the grain yield and N taken up by the biomass plotted in Figure 5. In general, the least N stress occurred in the year 2020 for N0, and low fertilization rates were observed, as the demand (uptake) by the plant was low, reflecting low productivity (grain yield). In contrast, in years with more favorable conditions, N stress increased slightly for low fertilization levels but all levelled out at a certain threshold (optimal fertilization) level.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Estimated average nitrogen stress index over the growing season (top row) and nitrogen use efficiency, agronomic N use efficiency ($AE_N$), partial factor productivity of applied N ($PFP_N$), and crop recovery efficiency of applied N ($RE_N$) with varying fertilization rates from 0 (N0) to 200 (N200) kg ha$^{-1}$.</figcaption></figure>

The calculated agronomic nitrogen efficiency coefficient ($AE_N$) is useful when investigating the economic portion of plants (grain yield). The overall trend of the $AE_N$ coefficient showed that 2021 had the highest values, indicating more favorable growing conditions compared to 2020 and 2022. Additionally, at higher N fertilization rates, $AE_N$ values increased up to a certain point before reaching a peak and then decreasing again. In 2020, $AE_N$ increased from 4.7 at a low fertilization of N10 to 14.3 at N100. Beyond N100, N efficiency started to decline, suggesting that additional N beyond N100 did not translate
to a proportional yield increase. In contrast, in 2021, AE<sub>N</sub> significantly increased from 7.4 at a low fertilization of N10 to 24.4 at N180, indicating a strong response of yield to N fertilization. In 2022, the calculated AE<sub>N</sub> coefficient increased from 4.1 at N10 to 13.1 at N170, followed by a gradual decline beyond this rate.

The partial factor productivity of applied nitrogen (PFP<sub>N</sub>) provides information on the ratio between produced barley grain yield and applied N. Therefore, it can indicate the crop's capability to convert N inputs into outputs (yield). The overall trends shown in Figure 6 indicate that the highest PFP<sub>N</sub> values were at the N10 rate across all years, suggesting low N rates were more efficient in terms of productivity per unit of applied N compared to higher rates. Additionally, all years showed almost the same pattern of PFP<sub>N</sub>, whereby for the year 2021, the low slope at medium and higher fertilization rates showed a higher PFP<sub>N</sub> compared to the years 2020 and 2022. This means that for the year 2021, more favorable climatic growth conditions led to a better N use efficiency, as indicated by the exceptionally high values compared to the other years. The low slope of PFP<sub>N</sub> at medium to higher fertilization rates suggests a diminishing return on N use efficiency at higher fertilization rates. Finally, the crop recovery efficiency of applied nitrogen (RE<sub>N</sub>) showed the proportion of applied N that was taken up and utilized by the crop. In general, the RE<sub>N</sub> trends were very similar to those observed for AE<sub>N</sub>, showing that in all three years, RE<sub>N</sub> values increased with higher N rates up to a peak and then began to decrease beyond at higher fertilization rates (see Figure 6). Similar to the AE<sub>N</sub> and PFP<sub>N</sub> coefficients, RE<sub>N</sub> indicated that the most favorable year for barley cultivation was 2021.

The highest RE<sub>N</sub> values in 2020 were calculated for the N110 rate, with RE<sub>N</sub> reaching 0.53. In 2021, RE<sub>N</sub> peaked at a much higher fertilization rate of N170 with a RE<sub>N</sub> of 0.88, while in 2022, it peaked at a fertilization rate of N170 with 0.49. However, all calculated N use efficiency indices did not account for the initial N reserves in the soil, which might confound the different years, even though care was taken to be independent on initial N stocks in the soil by adding a 5-year spin-up period in the simulation.

### 2.7.3. Estimated Water Stress and Water Use Efficiency

In order to evaluate the agronomic water use efficiency (WUE) for barley grown on the marginal sandy soil, the average water stress intensity for the entire barley growing season along with WUE was calculated and plotted in Figure 7 for the different N fertilization rates. In contrast to the N stress described above, the water stress intensity remained almost constant over the entire range of N fertilization rates and only differed slightly among the three years. For example, in 2020, water stress ranged from 0.75 to 0.77; in 2021, it was slightly lower with values ranging from 0.84 to 0.86; and in 2022, it showed the lowest values among all years at 0.73–0.75 over all fertilization rates. Water stress and its intensity mainly indicate the feedback between available water in the soil system replenished by precipitation and water losses by the plant due to transpiration. Classically, transpiration is coupled to photosynthetic activity or carbon dioxide assimilation for biomass growth as both carbon dioxide and water vapor are exchanged through the leaf stomata. Therefore, actual ET is an important indicator of crop water requirements and use [38]. However, for the sandy soils in the present study, a severe amount of water stress was given, no matter how much N fertilization was applied (see top panel, Figure 3).

The WUE values generally increased with higher N fertilization rates, then reached a peak, which corresponded to the maximum yield reached for the higher fertilization amounts and remained stable across all years (Figure 7). The highest WUE values (peaks) varied each year, reflecting the influence of yearly environmental conditions and the maximum yield reached in the individual years. In 2020, with the lowest total precipitation, WUE values for the N0 treatment were 2.88 kg ha<sup>-1</sup> mm<sup>-1</sup>, reaching a maximum
of 9.72 kg ha⁻¹ mm⁻¹ at the N120 fertilization rate. In 2021, with increased precipita-
tion, WUE values were higher, ranging from 2.98 to 16.40 kg ha⁻¹ mm⁻¹. Despite the
highest precipitation in 2022, WUE values were not the highest, ranging from 1.01 to
7.82 kg ha⁻¹ mm⁻¹. Accordingly, the total sum of precipitation over the growing season
did not solely determine the WUE rather than the distribution over the season. In 2022,
eight high rainfall events of 10.3, 11.3, 13.9, 14.5, 29.3, 30.8, 35.5, and 50.4 mm contributed
to a total of 196 mm of precipitation. Sandy soils, due to their large pores and high hy-
draulic conductivity at larger water contents, are known to exhibit also low water retention,
allowing water to infiltrate and drain quickly through the soil profile. Consequently, the in-
filtrated water may not be available for plant growth. As a result, a significant portion of the
2022 precipitation would be ineffective in maintaining barley growth and yield formation.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Estimated water availability index (α<sub>avg</sub>) (**top row**), water use efficiency (WUE) (**middle row**), and the correlation between the N availability index and WUE, with varying fertilization rates from 0 (N0) to 200 kg ha<sup>-1</sup> (N200), as well as N uptake and WUE (**bottom row**).</figcaption></figure>

Although in many cases, water stress can also lead to N stress, both stressors can
occur simultaneously under specific conditions, such as suboptimal water supply and
insufficient fertilization. Therefore, the relationships among WUE, N stress intensity, and N
uptake were analyzed by correlating WUE with the N stress intensity or plant N uptake
(see Figure 7). In all 3 years, a strong positive correlation was found between simulated
WUE and simulated N stress intensity (R² = 0.97–0.99), as well as between simulated WUE
and simulated N uptake into grains (R² = 0.99). Higher WUE values were associated with higher N availability, indicating that WUE increased when barley experienced less N stress during growth. A similar trend was observed between WUE and N uptake, emphasizing the interdependence of water and N efficiency in barley cultivation.

# 3. Materials and Methods

## 3.1. Site Description

According to the environmental stratification of Europe, Lithuania falls within the Nemoral climate zone. This climate zone is characterized by a continental and rather cool climate, as well as a relatively short vegetation period. The same climate zone encompasses the southern portion of Scandinavia, the Baltic states, and Belarus [39]. Lithuania's landscape territory has notable variations in terms of air temperature, precipitation pattern, and soil types. The mean annual air temperature ranges from 5.8 to 7.6 °C, while annual precipitation spans from 550 to 910 mm [40]. The predominant soil classification encompasses Luvisols, Cambisols, Gleysols, and Arenosols, constituting 28.5, 15.9, 14.6, and 13.2% of the land area, respectively.

The locations of the field trials fall into the agro-climatic zone IIIB of southeastern Lithuania [41], which has a pronounced seasonal variation with distinct spring, summer, autumn, and winter periods. This agro-climatic zone is relatively warm with moderate annual precipitation, but it is dominated by infertile (marginal) soils. The average annual air temperature in this zone is 7.2 °C and the total annual precipitation is 678 mm (mean value over the period 1991–2020). The main soil is Endocalcaric Eutric Brunic Arenosol (Geoabruptic, Aric, see Supplementary Figure S1), which is predominant in Lithuania [42]. The soil texture in the top soil horizon is sandy silt loam with an average soil organic carbon content of 1.34% (see Table 2). The main soil agrochemical characteristics were a slightly acidic pH, with high amounts of plant-available phosphorus (P₂O₅) and potassium (K₂O) levels (Table 2).

<table><caption>**Table 2.** Soil properties at the experimental site and for barley management information.</caption><thead><tr><th></th><th>2020</th><th>2021</th><th>2022</th></tr></thead><tbody><tr><td>Soil (FAO classification)</td><td colspan="3">Endocalcaric Eutric Brunic Arenosol (Geoabruptic, Aric)</td></tr><tr><td>Soil pH (1 N KCl extraction)</td><td>5.5</td><td>6.0</td><td>6.3</td></tr><tr><td>Soil P<sub>2</sub>O<sub>5</sub> (mg kg<sup>-1</sup>)<br>(Egner-Riehm-Domingo (A-L))</td><td>170</td><td>205</td><td>192</td></tr><tr><td>Soil K<sub>2</sub>O (mg kg<sup>-1</sup>)<br>(Egner-Riehm-Domingo (A-L))</td><td>324</td><td>174</td><td>180</td></tr><tr><td>Soil organic carbon (%)</td><td>1.34</td><td>1.16</td><td>1.21</td></tr><tr><td>Soil N total (%) (Kjeldahl)</td><td>0.103</td><td>0.101</td><td>0.103</td></tr><tr><td>Previous crop</td><td>Buckwheat</td><td>Barley</td><td>Barley</td></tr><tr><td>Barley cultivar</td><td>KWS Fantex</td><td>KWS Fantex</td><td>KWS Fantex</td></tr><tr><td>Barley seeding dates</td><td>27 April 2020</td><td>10 May 2021</td><td>04 May 2022</td></tr><tr><td>Seeding density (seeds per m<sup>2</sup>)</td><td>450</td><td>450</td><td>450</td></tr><tr><td>Plot size</td><td>3 m × 10 m = 30 m<sup>2</sup></td><td>3 m × 10 m = 30 m<sup>2</sup></td><td>3 m × 10 m = 30 m<sup>2</sup></td></tr><tr><td>Fertilization</td><td>27 April 2020—<br>N<sub>50</sub>P<sub>80</sub>K<sub>140</sub>;<br>8 June 2020—N<sub>50</sub></td><td>7 May 2021—N<sub>50</sub>P<sub>80</sub>K<sub>140</sub>;<br>11 June 2021—N<sub>50</sub></td><td>3 May 2022—N<sub>50</sub>P<sub>80</sub>K<sub>140</sub>;<br>6 June—N<sub>50</sub></td></tr><tr><td>Pesticides, fungicides, insecticides</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Barley harvesting</td><td>12 August 2020</td><td>13 August 2021</td><td>9 August 2022</td></tr></tbody></table>
## 3.2. Experimental Design

The field experiment was carried out at the Lithuanian Research Center for Agriculture and Forestry in the fields of the regional Vokė Branch (54.588029° N; 25.135752° W) in a 3 year period (2020–2022). The spring barley variety KWS Fantex was sown in late April/early May at a density of 450 m⁻² at 3–5 cm depth. The barley cultivar selected and used in the experiment was bred in Germany by the KWS Lochow GmbH Seed Company located in Bergen. It was chosen because of its high yield potential, high resistance to lodging, and very good resistance to net blotch, spot blotch, scald, and mildew, as well as double resistance to nematodes. Barley was annually harvested after the crops reached physiological maturity. The experimental design for the three years included various fertilization treatments. However, in this study, we focused on two treatments, (1) no fertilization (N0) and (2) N₁₀₀P₈₀K₁₄₀ (N100) plots with mineral fertilization, as these treatments provided the most contrasting nutrient conditions for this field experiment. Treatments were arranged in 4 randomized blocks with a plot size of 30 m², in which only 13.4 m² (1.68 × 8 m²) of the inner parts were harvested. The mineral fertilizers were in the form of ammonium nitrate (34.4-0-0), granular superphosphate (0-20-0), and potassium chloride (0-0-60). The full dose of P and K fertilizers was applied manually before sowing, while N application was split, i.e., N₅₀ was used before sowing and the N₅₀ rest was applied during the tillering growth stage of barley.

## 3.3. Plant Measurements

During the barley growth period, the development stages from emergence to maturity were assessed weekly using the BBCH scale [43]. The scale's first digit indicates the principal growth stages (e.g., 1 for leaf development and 2 for tillering), while the second digit provides details on a specific development step. A development stage for barley was designated when 50% or more of the plants across the plot had reached that particular stage.

From 2020 to 2022, the total aboveground biomass (TAB) of barley was periodically sampled throughout the growing season. Specifically, 6 TAB measurements were taken in 2020 at stages BBCH13, BBCH31, BBCH61, BBCH73, BBCH85, and at harvest. In 2021, 5 TAB measurements were conducted at BBCH21, BBCH32–37, BBCH75, BBCH85, and at harvest. In 2022, four TAB measurements occurred at BBCH23–25, BBCH47–59, BBCH75, and at harvest. For these measurements, samples were collected from each plot from a ground surface area of 0.5 m × 0.5 m (0.25 m²). The harvested plants were divided into their principal components (leaves, stems, and grains) as applicable on the date of sampling. Subsequently, the samples were dried in an oven at a temperature of 65 ± 5 °C until a constant weight was achieved to determine the dry biomass weights.

At barley physiological maturity (BBCH 92–99), the inner parts of experimental plots (1.68 m × 8 m = 13.4 m²) were harvested. Similar to the seasonal TAB measurements, at harvest, the plants were divided into their principal components and subjected to the same drying procedures. Barley leaves, stems, and grain subsamples taken at harvest were used to determine the N concentration in plant parts using the Kjeldahl method [44]. For each measured plot (N0 and N100) accumulated N uptake (based on dry matter) was calculated separately by multiplying the resulting nitrogen concentration (N%) by the corresponding yield by:

$$ N \text{ uptake} (\text{kg ha}^{-1}) = \frac{N\% \times \text{dry matter} (\text{kg ha}^{-1})}{100} \quad (1) $$

Additionally, three N use efficiency (NUE) coefficients were calculated: the agronomic N efficiency coefficient (AE_N), the partial factor productivity of applied N (PFP_N), and the
crop recovery efficiency of applied N (RE_N). These were calculated using the following equations, as described by Dobermann [45]:

$$AE_N = \frac{(Y_N - Y_0)}{F_N} \qquad (2)$$

$$PFP_N = Y_N - F_N \qquad (3)$$

$$REN = \frac{(U_N - U_0)}{F_N} \qquad (4)$$

where Y_N is the barley yield under the applied N rate (kg DM ha⁻¹), Y₀ is the barley yield (kg DM ha⁻¹) in a treatment with no N application, F_N is the amount of N applied (kg N ha⁻¹), U_N is the total barley N uptake in grain yield at maturity in a treatment with applied N fertilization (kg N ha⁻¹), and U₀ is the total barley N uptake in grain yield at maturity in a treatment without N fertilization (kg N ha⁻¹).

Leaf area was measured periodically using a simple HP printer scanner. Therefore, the plants were harvested from 0.25 m² area and separated into leaves, stem, and grains. In the next step, the leaves were scanned. The scanned images of the leaves were analyzed with the WinFOLIA (Regent Instruments Canada Inc., Québec City, QC, Canada) image analysis tool to calculate the leaf area index (LAI). Overall, four LAI measurements were conducted in 2020, and three each in 2021 and 2022.

## 3.4. Soil Measurements

### 3.4.1. Soil Chemical Properties

Before the start of each barley growing season, the nutrient status of the soil was evaluated. To conduct the soil analysis, composite samples were collected from the plowed layer (0–30 cm depth) at 12 distinct locations within the field. The chemical analyses included measurements of soil pH, soil organic carbon (SOC), total nitrogen (N_total), plant-available phosphorus (P₂O₅), potassium (K₂O), and bulk soil mineral nitrogen (SMN) as nitrate and ammonium. Additionally, two samples were taken from each of the five soil horizons for chemical analysis (Table 3). The pH of the soil samples was determined using an XS Instruments pH meter (Via della Meccanica, Italy) in a 1 M KCl solution with a 1:5 (vol/vol) ratio. The SOC content was measured through dry combustion using a Liqui TOC II instrument (Elementar, Langenselbold, Germany), with soil samples first being treated with HCl to eliminate the inorganic carbon fraction. These measurements were performed in triplicate for each sample. N_total was assessed using the Kjeldahl method with a semi-automatic Velp Scientifica™ UDK 139 instrument (VELP Scientifica, Usmate Velate, Italy), followed by manual titration of the sample with a 0.1 M NaOH solution. The soil P₂O₅ content was determined via a Shimadzu UV 1800 spectrophotometer, and the K₂O content was measured using flame emission spectroscopy with a JENWAY PFP7 flame photometer (Thermo Scientific, Cambridge, UK) at a 766 nm wavelength. Both nutrients were analyzed using the Egner-Riehm-Domingo (A-L) method. SMN was measured (ISO 14256-2:2005) using a spectrometric analyzer (Fiastar 5000, Foss Analytical AB, Höganäs, Sweden. All chemical analyses of the soil were performed at the Agrochemical Research Laboratory of the Lithuanian Research Centre for Agriculture and Forestry.
<table><caption>**Table 3.** Soil chemical properties determined before barley planting in 2020 at different soil horizons: Soil organic carbon (SOC), pH, total nitrogen (N<sub>total</sub>), plant-available phosphorus (P<sub>2</sub>O<sub>5</sub>), plant-available potassium (K<sub>2</sub>O), soil mineral nitrogen (SMN) including nitrate nitrogen (NO<sub>3</sub><sup>-</sup>-N), and ammonium nitrogen (NH<sub>4</sub><sup>+</sup>-N).</caption><thead><tr><th rowspan="2">Horizon<br>Description</th><th>SOC</th><th>pH</th><th>N<sub>total</sub></th><th>P<sub>2</sub>O<sub>5</sub></th><th>K<sub>2</sub>O</th><th>SMN<br>(NO<sub>3</sub><sup>-</sup>-N +<br>NH<sub>4</sub><sup>+</sup>-N)</th></tr><tr><th>%</th><th>-</th><th>%</th><th>mg<br>kg<sup>-1</sup></th><th>mg<br>kg<sup>-1</sup></th><th>kg ha<sup>-1</sup></th></tr></thead><tbody><tr><td>Ap (0–30 cm)</td><td>1.34</td><td>5.5</td><td>0.103</td><td>170</td><td>324</td><td>71.3 ± 11.5</td></tr><tr><td>B1 (30–50 cm)</td><td>0.46</td><td>5.8</td><td>0.013</td><td>120</td><td>153</td><td>–</td></tr><tr><td>B2 (50–78 cm)</td><td>0.54</td><td>6.4</td><td>0.014</td><td>259</td><td>122</td><td>–</td></tr><tr><td>2C(k) (78–105 cm)</td><td>0.37</td><td>7.6</td><td>0.003</td><td>195</td><td>43.5</td><td>–</td></tr><tr><td>2C (105–120 cm)</td><td>0.29</td><td>7.6</td><td>0.003</td><td>208</td><td>38</td><td>–</td></tr></tbody></table>

Horizon description: Ap = accumulative plaggic horizon in which decomposed organic material is being accumulated; B1 and B2 = illuvial horizons; 2C(k) and 2C = initial mineral horizons.

### 3.4.2. Soil Hydraulic Properties

In 2020, before the start of the barley growing season, a pit was manually dug in the central experimental location, and the soil was classified as an Endocalcaric Eutric Brunic Arenosol (Geoabruptic, Aric) with five main horizons according to WRB [42]. A photo of the profile is provided in Supplementary Figure S1. For soil physical characterization, 10 undisturbed soil samples of 250 cm³ were taken as duplicates from the horizons at 15–20, 40–45, 60–65, 90–95, and 110–115 cm. The soil hydraulic properties were determined by the evaporation method using the HYPROP® system (Meter Group, München, Germany) [46] in combination with the WP4® Dewpoint Potentiometer (Decagon Devices, WA, USA). Saturated hydraulic conductivity, K<sub>s</sub>, was measured using a falling head using the KSAT system (Meter Group, München, Germany). The estimated van Genuchten [47] soil hydraulic parameters for the soil profile up to a depth of 120 cm are listed in Table 4.

<table><caption>**Table 4.** Soil horizons and soil hydraulic parameters measured at the barley experimental location in 2020.</caption><thead><tr><th>Horizon<br>Description</th><th>PWP</th><th>θ<sub>r</sub><br>(cm<sup>3</sup> cm<sup>-3</sup>)</th><th>θ<sub>s</sub></th><th>α<br>(cm<sup>-1</sup>)</th><th>n<br>(-)</th><th>K<sub>s</sub><br>(cm day<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Ap (0–30 cm)</td><td>0.071</td><td>0.0145</td><td>0.419</td><td>0.011</td><td>1.395</td><td>58.87</td></tr><tr><td>B1 (30–50 cm)</td><td>0.043</td><td>0.0260</td><td>0.306</td><td>0.033</td><td>1.980</td><td>9.34</td></tr><tr><td>B2 (50–78 cm)</td><td>0.034</td><td>0.0205</td><td>0.285</td><td>0.040</td><td>3.647</td><td>5.88</td></tr><tr><td>2C(k) (78–105 cm)</td><td>0.056</td><td>0.0240</td><td>0.323</td><td>0.088</td><td>2.821</td><td>5.05</td></tr><tr><td>2C (105–120 cm)</td><td>0.026</td><td>0.0260</td><td>0.312</td><td>0.047</td><td>4.575</td><td>1.95</td></tr></tbody></table>

The soil hydraulic parameters according to van Genuchten [42] are: θ<sub>r</sub> = residual water content; θ<sub>s</sub> = saturated water content; a = inverse air entry pressure; n = shape parameter; and K<sub>s</sub> = saturated hydraulic conductivity estimated by calibration. PWP is the water content at the wilting point. Horizon description: Ap = accumulative plaggic horizon in which decomposed organic material is being accumulated; B1 and B2 = illuvial horizons; 2C(k) and 2C = initial mineral horizons.

Additionally, at all horizons, samples were taken to determine the soil texture and bulk density using the same undisturbed samples as used for the evaporation method (Table 5). Soil texture was analyzed according to DIN ISO 11277 [48] by wet sieving and the pipette method.
<table><caption>**Table 5.** Soil horizons, soil texture, and bulk density measured in 2020.</caption><thead><tr><th rowspan="2">Horizon<br>Description</th><th colspan="3">Particle Size %</th><th rowspan="2">Textural Class</th><th>Bulk<br>Density</th></tr><tr><th>Sand</th><th>Silt</th><th>Clay</th><th>(g cm<sup>-3</sup>)</th></tr></thead><tbody><tr><td>Ap (0–30 cm)</td><td>45.2</td><td>44.3</td><td>10.5</td><td>Sandy Silt Loam</td><td>1.47</td></tr><tr><td>B1 (30–50 cm)</td><td>88.0</td><td>7.4</td><td>4.6</td><td>Sand</td><td>1.59</td></tr><tr><td>B2 (50–78 cm)</td><td>81.9</td><td>11.4</td><td>6.7</td><td>Loamy Sand</td><td>1.64</td></tr><tr><td>2C(k) (78–105 cm)</td><td>90.9</td><td>5.9</td><td>3.2</td><td>Sand</td><td>1.76</td></tr><tr><td>2C (105–120 cm)</td><td>95.4</td><td>2.9</td><td>1.7</td><td>Sand</td><td>1.66</td></tr></tbody></table>

Note: Explanations of horizons are presented in Table 3, and sand was defined between 2000 and 63 µm according to WRB.

### 3.4.3. Soil Water Content

From 2020 to 2022, the soil volumetric water content (SWC) was monitored hourly at 5 different depths (15, 40, 60, 90, and 110 cm). For the SWC measurements, a wireless sensor network, SoilNet [49], was installed prior to seed bed preparation. In the AgroC model, daily averaged SWC values, aggregated from the hourly raw data, were used.

Barley water use efficiency (WUE) was calculated using the equation described by Howell [50]:

$$WUE = \frac{\text{Barley grain yield}}{\text{Water used to produce the yield}} \quad (5)$$

### 3.5. AgroC Model Setup and Calibration

For the simulation of barley growth, the AgroC model [21,22] was used. The AgroC model is a one-dimensional state-of-the-art agroecosystem model, which has been applied for different crops, such as winter wheat [22], maize [31,51], sugar beet [52], oat, barley [53], and hemp [20]. For the water flow, it solves the Richards equation, and the Mualem van Genuchten functions [47] were used for the soil hydraulic parameterization. More model details can be found in Herbst et al. [21] and Klosterhalfen et al. [22]. The detailed description of the AgroC N module is provided in Appendix A.

The AgroC model was set up to mimic the field experiments, in which the simulation domain was set to 150 cm depth and the 5 pedogenetic horizons, as listed in Table 3, were included. The upper boundary was defined by the atmospheric conditions with precipitation and given potential evapotranspiration (ET<sub>pot</sub>). The lower boundary was set to free drainage. Crop management (seeding, seeding density, harvest, and N-fertilization) followed the experimental data.

Model calibration was performed in two steps. First, phenology was calibrated using the approach suggested by Wallach et al. [54], in which the plant parameters, the input used [55] for barley simulation in western Germany, were taken as start parameters. Secondly, the soil hydraulics were optimized for each pedogenetic soil layer identified in the field using the measured soil water contents to be minimized in the objective function. To optimize the phenology and the soil hydraulic parameters, the Shuffled Complex Evolution (SCE-UA) algorithm [56,57] was used. As the soil hydraulic parameters were measured in the laboratory those parameters were used and only the saturated hydraulic conductivity K<sub>s</sub> was inversely estimated for each horizon, which minimized the number of parameters to assess. Please note that only one set of soil hydraulic parameters for the three plot locations, where the experiments were performed, has been optimized, as the small-scale variability in this area is typically low. Finally, the crop parameters were optimized manually to match observed biomasses, yields, and LAIs.

For all calibration steps, the first and third years (2020 and 2022) were selected, as those years are most contrasting in climate, fulfilling the requirements for calibration as
suggested by Wallach et al. [58]. The remaining second year (2021) was selected for model
validation. For the phenology and biomass (yield) calibration, all fertilization trial (N100)
data were used. To calibrate the nitrogen uptake and distribution within the plant organs,
both the N100 and no fertilization data for the N0 trials were used. As the soil nitrate and
ammonia concentrations were initialized by measured data, which had to be interpolated
and are known to be uncertain, the model was spun up for a 5-year period assuming the
same management as for the experimental years (same crop, same seeding and harvest
dates, and same fertilization). For the spin-up, the measured climatic data of the three years
were looped. An overview of the calibration steps and calibrated parameters is provided in
Table S1.

3.6. Climate Data

The weather data for this study were obtained from the Vilnius meteorological sta-
tion, located 5.3 km away from the barley experimental site. The daily data utilized for
simulations included precipitation (mm), maximum and minimum air temperature (Tmax
and Tmin) (°C), relative humidity (%), wind speed at 2 m height (m s⁻¹), and solar radia-
tion (J cm² day⁻¹). Using this data, the reference Penman–Monteith evapotranspiration
(ET₀) (mm day⁻¹) was calculated [59]. These weather data were sourced from the Lithua-
nian Hydrometeorological Service, through its data sharing service system (available at
https://www.meteo.lt, accessed on 8 January 2024).

3.7. AgroC Model Application

3.7.1. Prediction of Barley Yield Potential and Yield Gap

The impact assessment of the different environmental conditions on barley yield was
performed sequentially through simulations for the years 2020–2022. In a first step, the
potential yield was estimated by disabling the plant’s response mechanisms to nitrogen
and water, thereby removing abiotic stress factors and assuming ideal fertilization and
watering conditions. In a second step, the impact of water stress was simulated, in which
the nitrogen response was deactivated, isolating the effect of water only. In the next step,
the N impact on yields was assessed, by disabling the water stress and focusing solely on
the N impact. In the final step, actual barley yields were simulated by accounting for N and
water stress. The comparison of potential yields with those derived from the latter three
steps provided insight into the yield reductions caused by N stress, water stress, and the
combination of both stressors.

3.7.2. Evaluation of N Management Practice Scenarios

To evaluate the effect of increasing N fertilization rates from N0 (no fertilization) to N200 (200 kg N ha⁻¹) and climatic variables on barley total aboveground biomass, grain yield, water and N stress intensity, N uptake, and N leaching, the previously calibrated and validated AgroC model was used for simulations. The N fertilization scenarios were simulated for the same years in which the model was parameterized. For the different assumed N application levels of 10 to 200 kg ha⁻¹ per growing seasons, 19 simulations were performed with N application steps of 10 kg ha⁻¹ each (10, 20, 30, ..., 190, and 200 kg ha⁻¹). For fertilization, the same application dates as for the field experiments were used, and the assumed total N amounts were split in half for the first and second applications.

3.8. Statistical Analysis

The prediction capability of the model for both calibration and validation was tested by
the Willmott index of agreement (d), root mean squared error (RMSE), and BIAS, computed
as follows [60]:
$$d = 1 - \frac{\sum_{i=1}^{n} (Y_i - X_i)^2}{\sum_{i=1}^{n} (|Y_i - M0| + |X_i - M0|)^2} \quad (6)$$

$$RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - x_i)^2} \quad (7)$$

$$BIAS = 1/n \sum_{i=0}^{n} (y_i - x_i) \quad (8)$$

where $Y_i$ is the simulated value; $X_i$ is the measured value; MO is the average value of "n" measured values; and n is the number of measurements.

The index of agreement (d) developed by Willmott [60] as a standardized measure of the degree of model prediction error varies between 0 and 1. A value of 1 indicates a perfect match between observation and the model and 0 indicates no agreement at all. The BIAS measures the average difference between measured and simulated values. A positive BIAS indicates model under-prediction, and a negative BIAS indicates over-prediction. A lower RMSE indicates a better fit to the observed data compared to larger values.

## 4. Conclusions

In this study, the AgroC model was calibrated and validated using a multiyear barley field experiment on marginal sandy soil under varying fertilization and weather conditions. After calibration, the model accurately simulated barley development, grain yield, and N uptake. It was then used to predict N leaching, assess N and water stress effects on yields, and analyze potential yield and yield gaps beyond what experimental data alone could provide.

This study provided experimental and modeling evidence of potential barley grain yield levels of 4.8–6.02 t ha⁻¹ (dry weight) on marginal soil at a latitude of 54° N in the Baltic region. However, a significant barley yield gap was identified, primarily due to limited water and nitrogen availability in these poor sandy soils. To address this, farmers should integrate fertilization and irrigation management practices. The combined effects of water and full N stress on the sandy soil reduced barley yields by more than 80% for all three experimental years pointing to the large stress pressure at this site. To explore optimal fertilization rates, a systemic study was performed based on the calibrated model, in which the modeling results for different fertilization scenarios suggested that higher N fertilization has the potential to increase actual yields, although this increase varied between the years making general recommendations complicated. The highest barley yields were achieved with N rates between 110 and 180 kg ha⁻¹. However, with increasing N fertilization rates, fertilizer efficiency declined, while nitrate leaching intensified. Therefore, the optimal and economically beneficial N fertilization rate for barley on this sandy soil should be between 100 and 120 kg ha⁻¹.

This study contributed to the ongoing efforts to develop effective strategies for reducing yield gaps, ensuring food security, and promoting sustainable agricultural development in the Baltic–Nordic region.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/plants14050704/s1, Figure S1: Soil profile and its main horizons; Table S1: Overview of the calibration steps and the calibrated parameters for the AgroC model.

**Author Contributions:** Conceptualization, R.Ž., L.W., M.H. methodology, R.Ž., L.W., M.H.; software, R.C., L.W., M.H.; formal analysis, R.Ž., L.W., M.H.; investigation, R.Ž., V.P., E.B.; writing—original draft preparation, R.Ž., L.W.; writing—review and editing, R.C., L.W., M.H., E.L.-H., W.S., P.S., V.P., M.M., F.R., E.B., J.V., B.R., A.P.; visualization, R.Ž. All authors have read and agreed to the published version of the manuscript.
**Funding:** The barley experimental data were collected as part of the BioFoodonMars project, which was conducted from 2020 to 2023. This project was funded by the EU-FACCE-SURPLUS initiative within the framework of FACCE-JPI. The modeling portion of the study was carried out during a postdoctoral fellowship. Additionally, this project received funding from the Research Council of Lithuania (LMTLT), under agreement No S-PD-24-64.

**Data Availability Statement:** The original contributions presented in this study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

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

## Appendix A

The AgroC nitrogen module accounts for the carbon turnover-linked mineralization and immobilization of soil organic N, nitrification, denitrification, urea hydrolysis, ammonia volatilization, N₂O emission, N fixation, and crop demand-driven uptake by roots. In conjunction with the solute transport routine, ammonium sorption, convective-dispersive transport, and leaching of ammonium and nitrate can be accounted for (see Figure A1).

<figure><img src="image_9.png" /><figcaption>Figure A1. AgroC nitrogen cycle as a combination of RothC [61] and SWATNIT [62].</figcaption></figure>

### Appendix A.1. Mineralization of Soil Organic Nitrogen

The potential rate constants of the organic matter pools are specified as inputs in the carbon turnover section. The actual turnover rates are consistently used for the decomposition of organic C, N, and P. Actual turnover rates are calculated from the potential rate constants independent of soil temperature, water content (or pressure head), and aeration. The organic nitrogen pool structure follows the structure of the RothC pool concept applied for the organic C turnover, except for the BIO pool. This pool explicitly accounts for the carbon stored in the biomass and is for the organic nitrogen turnover replaced with a fixed biomass C/N ratio provided as model input. The decomposition of organic N in the decomposable (DPM) and resistant plant material (RPM) as well as in the humus pool (HUM) is estimated as follows:

$$ \frac{\partial N_{DPM}}{\partial t} = \left[ (1 - f_h) f_e \frac{1}{r_o} - \frac{N_{DPM}}{C_{DPM}} \right] C_{DPM} \lambda_{DPM} \quad (A1) $$
where N<sub>DPM</sub> [M N L<sup>-3</sup> soil] is the nitrogen content in the decomposable plant material DPM pool, C<sub>DPM</sub> [M C L<sup>-3</sup> soil] is the carbon content of the DPM pool, λ<sub>DPM</sub> [T<sup>-1</sup>] is the DPM turnover rate, and r<sub>o</sub> [-] represents the biomass C/N ratio. For consistency with the RothC estimated carbon turnover, the dimensionless turnover efficiency f<sub>e</sub> [-] is computed as 1-b/(b + 1), and the dimensionless humification factor f<sub>h</sub> is equal to 0.54 (=RothC humification factor). A similar nitrogen turnover loop is assumed for the resistant plant material pool:

$$ \frac{\partial N_{RPM}}{\partial t} = \left[ (1 - f_h) f_e \frac{1}{r_o} - \frac{N_{RPM}}{C_{RPM}} \right] C_{RPM} \lambda_{RPM} \quad (A2) $$

where N<sub>RPM</sub> [M N L<sup>-3</sup> soil] is the nitrogen content in the RPM pool, C<sub>RPM</sub> [M C L<sup>-3</sup> soil] is the carbon content of the RPM pool, and λ<sub>RPM</sub> [T<sup>-1</sup>] is the RPM turnover rate. The net increase rate of nitrogen in the humus pool N<sub>HUM</sub> [M N L<sup>-3</sup> soil] is given by:

$$ \frac{\partial N_{HUM}}{\partial t} = \frac{f_e f_h}{r_o} (C_{DPM} \lambda_{DPM} + C_{RPM} \lambda_{RPM}) - C_{HUM} \lambda_{HUM} \quad (A3) $$

where C<sub>HUM</sub> [M C L<sup>-3</sup> soil] is the carbon content of the HUM pool and λ<sub>HUM</sub> [T<sup>-1</sup>] is the HUM turnover rate. The C and N turnover may result in net production or consumption of ammonium, which depends on the biomass C/N ratio and the C/N ratio of the three pools organic matter pools. The decrease in soil organic N in the three pools due to mineralization equals the increase in nitrogen in the form of ammonium N<sub>NH4</sub> [M N L<sup>-3</sup> soil]:

$$ \frac{\partial N_{NH4}}{\partial t} = \left[ \frac{N_{DPM}}{C_{DPM}} - \frac{f_e}{r_o} \right] C_{DPM} \lambda_{DPM} + \left[ \frac{N_{RPM}}{C_{RPM}} - \frac{f_e}{r_o} \right] C_{RPM} \lambda_{RPM} + C_{HUM} \lambda_{HUM} \quad (A4) $$

Ammonium-N is converted to the mass of NH<sub>4</sub> and to a liquid phase concentration to be used as a source/sink term [M NH<sub>4</sub> L<sup>-3</sup> water T<sup>-1</sup>] in the solute transport equation. When the right side of the equation above is positive, mineralization occurs. In case the right side is negative, immobilization occurs, and nitrogen mass is transferred from the liquid phase to the organic matter nitrogen pools. If immobilization occurs and the liquid phase concentration of NH<sub>4</sub> is not sufficient to fulfill the mineralization N demand, liquid phase NO<sub>3</sub> is transferred to the organic N pools. If the NO<sub>3</sub> supply to the mineralization process is not sufficient, λ<sub>DPM</sub> and λ<sub>RPM</sub> are reduced by the supply/demand ratio. For the supply from the liquid phase N, the N<sub>DPM</sub> pool is prioritized to the N<sub>RPM</sub> pool, i.e., the remaining liquid phase ammonium and/or nitrate may be consumed by the N<sub>DPM</sub> decomposition first. The potentially decreased λ<sub>DPM</sub> and λ<sub>RPM</sub> rates are also used in the organic C and P mineralization module. Thus, a nitrogen gap in the organic N mineralization of one or both plant material pools will likewise reduce the decomposition of the respective plant material C and P pools. Equations (A1)–(A4) are solved numerically by application of the Euler integration method.

## Appendix A.2. Nitrification, Ammonia Volatilization and Urea Hydrolysis

In the nitrification process, ammonium is transformed to nitrate by heterotrophic microorganisms. A first order kinetic equation, based on an actual nitrification rate is applied, in which the actual nitrification rate is computed from a potential nitrification rate λ<sub>nitri</sub> [T<sup>-1</sup>] specified as input. The potential nitrification rate is scaled with the rate modifiers accounting for soil temperature and water content also used for the organic decomposition rate scaling. Ammonium sorbed to the soil matrix is not transformed to nitrate; this only occurs for the dissolved ammonium.

Ammonia volatilization is again described by a first order kinetic approach, in which the actual volatilization rate is computed from the potential volatilization rate, λ<sub>vol</sub> [T<sup>-1</sup>], scaled by the rate modifiers for soil temperature and water content.
Hydrolysis of urea, as a basic constituent in many mineral N fertilizers, is also described with a first order kinetic reaction. Again, the soil water content and soil temperature affect the potential hydrolysis rate, $\lambda_{hyd}$ [T⁻¹].

## Appendix A.3. Denitrification

Denitrification is simulated as a first order kinetic reaction, in which the liquid phase nitrate concentration $NO_3$ [M NO₃ L⁻³ water] is denitrified according to the actual denitrification rate $\lambda_{a,denit}$ [T⁻¹]:

$$ \frac{\partial NO_3}{\partial t} = NO_3\lambda_{a,denit} \quad (A5) $$

The actual denitrification rate is computed from the potential denitrification rate, $\lambda_{denit}$ [T⁻¹], given as input, and it is scaled with the common rate modifier for soil temperature. However, the rate modifier for soil water content, $f_{w,denit}$ [-], differs from the soil water content rate modifier used for all the other turnover rates in order to account for the effect of soil aeration:

$$ \begin{gathered} f_{w,denit} = 0 \quad \text{for } \Theta_e \le \theta_d \\ f_{w,denit} = \left( \frac{\Theta_e - 0.8}{0.2} \right)^2 \quad \text{for } \Theta_e > \theta_d \end{gathered} \quad (A6) $$

where $\theta_d$ [L³ L⁻³] is a threshold water content, set to 0.8 according to Aulakh et al. [63], $\Theta_e$ is the effective saturation [-], and $\theta_s$ [L³ L⁻³] is the water content at saturation.

## Appendix A.4. Root Uptake of Nitrogen

The uptake scheme follows the approach [64], which was modified [65].

The potential uptake (=demand), $N_p$ [M N L⁻²], is computed from a dry matter nitrogen concentration for each crop organ, as provided in the plant parameter input file against development stage DVS. Total potential uptake rate $dN_p/dt$ [M N L⁻² T⁻¹] is defined as the sum of the following organ-specific potential uptake rates:

$$ \frac{\partial N_p}{\partial t} = \frac{\partial N_{lea}}{\partial t} + \frac{\partial N_{ste}}{\partial t} + \frac{\partial N_{roo}}{\partial t} + \frac{\partial N_{crn}}{\partial t} \quad (A7) $$

where $N_{lea}$, $N_{ste}$, $N_{roo}$, and $N_{crn}$ [M N L⁻²] represent the nitrogen demand of leaves, stems, storage organs, roots, and crowns, respectively. The potential uptake rate of each plant organ is computed as:

$$ \frac{\partial N_i}{\partial t} = W_i * XNC_i - N_{act,i} \quad (A8) $$

where $W_i$ [M DM L⁻²] is the accumulated dry matter, $XNC_i$ [M N M⁻¹ DM] is the potential dry matter nitrogen concentration of the crop organs, $N_{act,i}$ [M N L⁻²] is the nitrogen accumulated in the organ, and index i loops over leaves, stems, storage organs, and roots. Uptake is supposed to be convective, with the root water uptake diffusive. The potential convective uptake rate [M N T⁻¹] is given by:

$$ \frac{\partial N_{conv}}{\partial t} = \int_{0}^{R_d} S_w * C_m * dx \quad (A9) $$

where $R_d$ is the actual root depth [L], $dx$ is the element thickness [L], $C_m$ is the liquid phase concentration of either nitrate [M NO₃ L⁻³ water] or ammonium [M NH₄ L⁻³ water], and $S_w$ [L³ L⁻³ T⁻¹] is the root water uptake sink term. The nitrogen uptake is simulated as a sink term in the solute transport routine, in which the uptake could never exceed the available liquid phase nitrogen mass. Ammonium will only be taken up when the crop nitrogen demand exceeds the nitrate supply. Only when convective uptake is smaller than
the potential uptake is the potential diffusive uptake calculated as the difference between potential and convective uptake. Actual diffusive uptake N<sub>dif</sub> [M N L<sup>-2</sup>] is calculated as:

$$ \frac{\partial N_{dif}}{\partial t} = \int_{0}^{R_d} \frac{2 * \pi * R_{dens} * RORAD * \tau * C_m * \theta}{D_0} dx \quad (A10) $$

where R<sub>dens</sub> [L L<sup>-3</sup>] is the root density at a specific depth, RORAD [L] is the mean root radius, τ [L<sup>2</sup> T<sup>-1</sup>] is the solute diffusion coefficient for a given water content, D<sub>0</sub> [L<sup>-1</sup>] is the travel distance resistance between the bulk soil solution and roots, C<sub>m</sub> is the liquid phase concentration of either nitrate [M NO<sub>3</sub> L<sup>-3</sup> water] or ammonia [M NH<sub>4</sub> L<sup>-3</sup> water], and dx is the depth increment [L]. Please note that the uptake rates N<sub>dif</sub> and N<sub>conv</sub> are related to soil volume and not to water volume. Root density over depth is calculated by scaling the root density at the surface W<sub>0dens</sub> [L L<sup>-3</sup>] with the relative root density RRD [-].

The organ-specific actual root nitrogen uptake, dRNUP<sub>act,i</sub>/dt [M N L<sup>-2</sup>], in which index i again loops over leaves, stems, storage organs, and roots is given as:

$$ \frac{\partial RNUP_{act,lea}}{\partial t} = \frac{\partial N_{lea}}{\partial t} * \frac{\frac{\partial N_{conv} + \partial N_{dif}}{\partial t}}{\frac{\partial N_p}{\partial t}} - \frac{\partial N_{so}}{\partial t} * \left( \frac{W_{lea}}{W_{lea} + W_{ste}} \right) \quad (A11) $$

$$ \frac{\partial RNUP_{act,ste}}{\partial t} = \frac{\partial N_{ste}}{\partial t} * \frac{\frac{\partial N_{conv} + \partial N_{dif}}{\partial t}}{\frac{\partial N_p}{\partial t}} - \frac{\partial N_{so}}{\partial t} * \left( \frac{W_{ste}}{W_{lea} + W_{ste}} \right) \quad (A12) $$

$$ \frac{\partial RNUP_{act,roo}}{\partial t} = \frac{\partial N_{roo}}{\partial t} * \frac{\frac{\partial N_{conv} + \partial N_{dif}}{\partial t}}{\frac{\partial N_p}{\partial t}} \quad (A13) $$

$$ \frac{\partial RNUP_{act,so}}{\partial t} = \frac{\partial N_{so}}{\partial t} = (W_{so} * XNC_{so} - N_{so}) * FN_{def} \quad (A14) $$

$$ FN_{def} = 1 - \sqrt{1 - N_{RED}^2} \quad (A15) $$

where FN<sub>def</sub> [-] describes the reduction of N transfer to the storage organs according to a potential limitation of N uptake.

## Appendix A.5. N Stress

The overall reduction of growth, N<sub>RED</sub> [-], according to nitrogen limitation, is described with a dimensionless reduction factor, varying between 0 and 1:

$$ N_{RED} = \left( \frac{ANCL - RLNCL}{RMNCL - RLNCL} \right) \quad (A16) $$

where ANCL [kg N kg<sup>-1</sup> DM] is the actual tissue nitrogen concentration of the leaves, RMNCL [kg N kg<sup>-1</sup> DM] is a leaf nitrogen threshold value for unrestricted growth assumed to be 50% of XNC<sub>LE</sub>, and RLNCL represents the leaf nitrogen threshold concentration below which growth ceases. The latter is fixed to 0.005 [M N M<sup>-1</sup> DM].

## Appendix A.6. Water Uptake Stress

Potential transpiration is distributed over the soil depth according to the relative root density distribution to provide the potential sink term of root water uptake over soil depth. The depth-specific actual root water uptake is computed by scaling the potential root water
uptake with the reduction factor α (-) independent of soil pressure head h ([L]) following the approach [66]:

$$ \alpha(h) = \begin{cases} \frac{h_0-h}{h_0-h_1} & h_0 \le h \le h_1 \\ 1 & \text{for } h_1 \le h \le h_2 \\ 10^{-\frac{h_2-h}{h_3}} & h_2 \le h \le h_3 \end{cases} \qquad (A17) $$

where h₀, h₁, h₂, and h₃ ([L]) are prescribed threshold pressure heads [67], which are plant dependent. Integrating the actual root water uptake over depth provides the actual transpiration Tₐ ([L T⁻¹]). The reduction of stomatal conductance due to water stress corresponds to the ratio between actual and potential transpiration Tₐ/Tₚ and is referred to as aₐₘg.

## Appendix A.7. Effect of Stress on Growth

The effect of stressors, each varying between 0 and 1, on photosynthesis is assumed to follow the minimum concept:

$$ GA_{t,a} = GA_{t,p} * \text{Min}(\alpha_{avg}, N_{RED}, P_{RED}) \qquad (A18) $$

where GAₜₐ [M CH₂O L⁻² T⁻¹] is the actual total gross assimilation, GAₜₚ [M CH₂O L⁻² T⁻¹] is the potential total gross assimilation, aₐₘg [-] is the stress factor accounting for water stress, and PₐₘED [-] is the stressor potentially accounting for limitations in phosphorus uptake.

## References

1. IPCC. Climate Change 2014: Impacts, Adaptation, and Vulnerability; Intergovernmental Panel on Climate Change, IPCC Secretariat: Geneva, Switzerland, 2014.

2. United Nations Department of Economic and Social Affairs. World Population Prospects 2024: Summary of Results. Available online: https://desapublications.un.org/publications/world-population-prospects-2024-summary-results (accessed on 10 September 2024).

3. Schröder, P.; Mench, M.; Povilaitis, V.; Rineau, F.; Rutkowska, B.; Schloter, M.; Szulc, W.; Žydelis, R.; Loit, E. Relaunch cropping on marginal soils by incorporating amendments and beneficial trace elements in an interdisciplinary approach. *Sci. Total Environ.* **2022**, *803*, 149844. [CrossRef] [PubMed]

4. Appiah, M.; Abdulai, I.; Schulman, A.; Moshelion, M.; Dewi, E.S.; Daszkowska-Golec, A.; Bracho-Mujica, G.; Rötter, R.P. Drought response of water-conserving and non-conserving spring barley cultivars. *Front. Plant Sci.* **2023**, *14*, 1247853. [CrossRef] [PubMed]

5. Wang, E.; Martre, P.; Zhao, Z.; Ewert, F.; Maiorano, A.; Rötter, R.P.; Kimball, B.A.; Ottman, M.J.; Wall, G.W.; White, J.W.; et al. The uncertainty of crop yield projections is reduced by improved temperature response functions. *Nat. Plants* **2017**, *3*, 17102. [CrossRef] [PubMed]

6. Wang, E.; Prentice, I.C.; Keenan, T.F.; Davis, T.W.; Wright, I.J.; Cornwell, W.K.; Evans, B.J.; Peng, C. Towards a universal model for carbon dioxide uptake by plants. *Nat. Plants* **2017**, *3*, 734–741. [CrossRef] [PubMed]

7. Vogel, E.; Donat, M.G.; Alexander, L.V.; Meinshausen, M.; Ray, D.K.; Karoly, D.; Meinshausen, N.; Frieler, K. The effects of climate extremes on global agricultural yields. *Environ. Res. Lett.* **2019**, *14*, 054010. [CrossRef]

8. Žydelis, R.; Dechmi, F.; Isla, R.; Weihermüller, L.; Lazauskas, S. CERES-Maize model performance under mineral and organic fertilization in nemoral climate conditions. *Agron. J.* **2021**, *113*, 2474–2490. [CrossRef]

9. Malik, W.; Isla, R.; Dechmi, F. DSSAT-CERES-maize modelling to improve irrigation and nitrogen management practices under Mediterranean conditions. *Agric. Water Manag.* **2019**, *213*, 298–308. [CrossRef]

10. Povilaitis, A.; Šileika, A.; Deelstra, J.; Gaigalis, K.; Baigys, G. Nitrogen losses from small agricultural catchments in Lithuania. *Agric. Ecosyst. Environ.* **2014**, *195*, 54–64. [CrossRef]

11. Bijay-Singh; Craswell, E. Fertilizers and nitrate pollution of surface and ground water: An increasingly pervasive global problem. *SN Appl. Sci.* **2021**, *3*, 518. [CrossRef]

12. Yawson, D.O.; Adu, M.O.; Armah, F.A. Impacts of climate change and mitigation policies on malt barley supplies and associated virtual water flows in the UK. *Sci. Rep.* **2020**, *10*, 376. [CrossRef] [PubMed]
13. Dawson, I.K.; Russell, J.; Powell, W.; Steffenson, B.; Thomas, W.T.B.; Waugh, R. Barley: A translational model for adaptation to climate change. *New Phytol.* **2015**, *206*, 913–931. [CrossRef] [PubMed]

14. FAOSTAT. *Statictical Database*; Food and Agriculture Organization of the United Nations: Rome, Italy, 2023. Available online: https://www.fao.org/faostat/en/#data/QCL (accessed on 12 February 2025).

15. Olesen, J.E.; Trnka, M.; Kersebaum, K.C.; Skjelvåg, A.O.; Seguin, B.; Peltonen-Sainio, P.; Rossi, F.; Kozyra, J.; Micale, F. Impacts and adaptation of European crop production systems to climate change. *Eur. J. Agron.* **2011**, *34*, 96–112. [CrossRef]

16. Khanna, M.; Chen, L.; Basso, B.; Cai, X.; Field, J.L.; Guan, K.; Jiang, C.; Lark, T.L.; Richard, T.L.; Spawn-Lee, S.A.; et al. Redefining marginal land for bioenergy crop production. *GCB Bioenergy* **2021**, *13*, 1590–1609. [CrossRef]

17. Gerwin, W.; Repmann, F.; Galatsidas, S.; Vlachaki, D.; Gounaris, N.; Baumgarten, W.; Volkmann, C.; Keramitzis, D.; Kiourtsis, F.; Freese, D. Assessment and quantification of marginal lands for biomass production in Europe using soil-quality indicators. *Soil* **2018**, *4*, 267–290. [CrossRef]

18. Persson, T.; Höglind, M.; Wallsten, J.; Nadeau, E.; Huang, X.; Rustas, B. Adapting the grassland model BASGRA to simulate yield and nutritive value of whole-crop barley. *Eur. J. Agron.* **2024**, *154*, 127075. [CrossRef]

19. Kersebaum, K.C.; Boote, K.J.; Jorgenson, J.S.; Nendel, C.; Bindi, M.; Frühauf, C.; Gaiser, T.; Hoogenboom, G.; Kollas, C.; Olesen, J.E.; et al. Analysis and classification of data sets for calibration and validation of agro-ecosystem models. *Environ. Model. Softw.* **2015**, *72*, 402–417. [CrossRef]

20. Žydelis, R.; Herbst, M.; Weihermüller, L.; Ruzgas, R.; Volungevičius, J.; Barčauskaitė, K.; Tilvikienė, V. Yield potential and factor influencing yield gap in industrial hemp cultivation under nemoral climate conditions. *Eur. J. Agron.* **2022**, *139*, 126576. [CrossRef]

21. Herbst, M.; Hellebrand, H.J.; Bauer, J.; Huisman, J.A.; Šimůnek, J.; Weihermüller, L.; Graf, A.; Vanderborght, J.; Vereecken, H. Multiyear heterotrophic soil respiration: Evaluation of a coupled CO₂ transport and carbon turnover model. *Ecol. Model.* **2008**, *214*, 271–281. [CrossRef]

22. Klosterhalfen, A.; Herbst, M.; Weihermüller, L.; Graf, A.; Schmidt, M.; Stadler, A.; Schneider, K.; Subke, J.A.; Huisman, J.A.; Vereecken, H. Multi-site calibration and validation of a net ecosystem carbon exchange model for croplands. *Ecol. Model.* **2017**, *363*, 137–156. [CrossRef]

23. Official Statistics Portal. Agricultural Production Methods and Irrigation. 2024. Available online: https://osp.stat.gov.lt/zus2020-rezultatai/zemes-ukio-gamybos-metodai-ir-drekinimas (accessed on 12 February 2025).

24. Wang, N.; Wang, J.; Wang, E.; Yu, Q.; Shi, Y.; He, D. Increased uncertainty in simulated maize phenology with more frequent supra-optimal temperature under climate warming. *Eur. J. Agron.* **2015**, *71*, 19–33. [CrossRef]

25. Usowicz, B.; Lipiec, J. Spatial variability of soil properties and cereal yield in a cultivated field on sandy soil. *Soil Tillage Res.* **2017**, *174*, 241–250. [CrossRef]

26. Castaldi, F.; Koparan, M.H.; Wetterlind, J.; Žydelis, R.; Vinci, I.; Savas, A.O.; Kivrak, C.; Tuncay, T.; Volungevičius, J.; Obber, S.; et al. Assessing the capability of Sentinel-2 time-series to estimate soil organic carbon and clay content at local scale in croplands. *ISPRS J. Photogramm. Remote Sens.* **2023**, *199*, 40–60. [CrossRef]

27. Wiegmann, M.; Maurer, A.; Pham, A.; March, J.T.; Al-Abdallat, A.; Thomas, W.T.B.; Bull, H.J.; Shahid, M.; Eglinton, J.; Baum, M.; et al. Barley yield formation under abiotic stress depends on the interplay between flowering time genes and environmental cues. *Sci. Rep.* **2019**, *9*, 6397. [CrossRef] [PubMed]

28. Baik, B.K.; Newman, C.W.; Newman, R.K. *Barley: Production, Improvement, and Uses*; Ullrich, S.E., Ed.; Wiley-Blackwell: Oxford, UK, 2011. [CrossRef]

29. Schils, R.; Olesen, J.E.; Kersebaum, K.C.; Rijk, B.; Oberforster, M.; Kalyada, V.; Khitrykau, M.; Gobin, A.; Kirchev, H.; Manolova, V.; et al. Cereal yield gaps across Europe. *Eur. J. Agron.* **2018**, *101*, 109–120. [CrossRef]

30. Hoffman, M.P.; Haakana, M.; Asseng, S.; Höhn, J.G.; Palosuo, T.; Ruiz-Ramos, M.; Fronzek, S.; Ewert, F.; Gaiser, T.; Kassie, B.T.; et al. How does inter-annual variability of attainable yield affect the magnitude of yield gaps for wheat and maize? An analysis at ten sites. *Agric. Syst.* **2018**, *159*, 199–208. [CrossRef]

31. Žydelis, R.; Weihermüller, L.; Herbst, M.; Klosterhalfen, A.; Lazauskas, S. A model study on the effect of water and cold stress on maize development under nemoral climate. *Agric. For. Meteorol.* **2018**, *263*, 169–179. [CrossRef]

32. Van Grinsvel, H.J.M.; Ebanyat, P.; Glendining, M.; Gu, B.; Hijbeek, R.; Lam, S.K.; Lassaletta, L.; Mueller, N.D.; Pacheco, F.S.; Quemada, M.; et al. Establishing long-term nitrogen response of global cereals to assess sustainable fertilizer rates. *Nat. Food* **2022**, *3*, 122–132. [CrossRef] [PubMed]

33. EC. *Council Directive 91/676/EEC of 12 December 1991 Concerning the Protection of Waters Against Pollution Caused by Nitrates from Agricultural Sources*; European Commission (EC): Brussels, Belgium, 1991.

34. Valkama, E.; Salo, T.; Esala, M.; Turtola, E. Nitrogen balances and yields of spring cereals as affected by but nitrogen fertilization in northern conditions: A meta-analysis. *Agric. Ecosyst. Environ.* **2013**, *164*, 1–13. [CrossRef]

35. *Tręšimo Plano Sudarymas*; Lietuvos Žemės Ūkio Konsultavimo Tarnyba: Kėdainiu, Lithuania, 2002; ISBN 9986-540-83-6.

36. World Health Organization. *Guidelines for Drinking-Water Quality*; WHO: Geneva, Switzerland, 2012.
37. Abascal, E.; Gómez-Coma, L.; Ortiz, I.; Ortiz, A. Global diagnosis of nitrate pollution in groundwater and review of removal technologies. *Sci. Total Environ.* **2022**, *810*, 152233. [CrossRef]

38. Goyal, M.R.; Harmsen, E.W. *Evapotranspiration*; Apple Academic Press: New York, NY, USA, 2013. [CrossRef]

39. Metzger, M.J.; Shkaruba, A.D.; Jongman, R.H.G.; Bunce, R.G.H. *Description of the European Environmental Zones and Strata*; Alterra Report 2281; Alterra: Wageningen, The Nederland, 2012.

40. Kavaliauskas, A.; Žydelis, R.; Castaldi, F.; Auškalnienė, O.; Povilaitis, V. Predicting Maize Theoretical Methane Yield in Combination with Ground and UAV Remote Data Using Machine Learning. *Plants* **2023**, *12*, 1823. [CrossRef]

41. Bukantis, A. Agroclimatic Zoning. In *Lithuanian National Atlas*; National Land Service under the Ministry of Agriculture: Vilniaus, Lithuania, 2009.

42. WRB. *World Reference Base for Soil Resources 2014*. *International Soil Classification System for Naming Soils and Creating Legends for Soil Maps*; World Soil Resources Reports No. 106; FAO: Rome, Italy, 2014.

43. Meier, U. Growth Stages of Mono- and Dicotyledonous Plants. BBCH Monograph. Julius Kühn-Institut (JKI) Quedlinburg. 2018. Available online: https://www.julius-kuehn.de/media/Veroeffentlichungen/bbch%20epaper%20en/page.pdf (accessed on 10 January 2023).

44. Mažeika, R.; Arbačiauskas, J.; Masevičienė, A.; Narutytė, I.; Šumskis, D.; Žičkienė, L.; Rainys, K.; Drapanauskaitė, D.; Staugaitis, G.; Baltrusaitis, J. Nutrient Dynamics and Plant Response in Soil to Organic Chicken Manure-Based Fertilizers. *Waste Biomass Valorization* **2020**, *12*, 371–382. [CrossRef]

45. Dobermann, A.R. Nitrogen use efficiency-state of the art. In *Proceedings of the IFA International Workshop on Enhanced-Efficiency Fertilizers*, Frankfurt, Germany, 28–30 June 2005.

46. Schindler, U.; Durner, W.; von Unold, G.; Mueller, L. Evaporation method for measuring unsaturated hydraulic properties of soils: Extending the measurement range. *Soil Sci. Soc. Am. J.* **2010**, *74*, 1071–1083. [CrossRef]

47. van Genuchten, M.T. A close-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Am. J.* **1980**, *44*, 892–898. [CrossRef]

48. Müller, H.W.; Dohrmann, R.; Klosa, D.; Rehder, S.; Eckelmann, W. Comparison of two procedures for particle-size analysis: Köhn pipette and X-ray granulometry. *J. Plant Nutr. Soil Sci.* **2009**, *172*, 172–179. [CrossRef]

49. Bogena, H.R.; Huisman, J.A.; Schilling, B.; Weuthen, A.; Vereecken, H. Effective Calibration of Low-Cost Soil Water Content Sensors. *Sensors* **2017**, *17*, 208. [CrossRef]

50. Howell, T.A. Enhancing Water Use Efficiency in Irrigated Agriculture. *Agron. J.* **2001**, *93*, 281–289. [CrossRef]

51. Žydelis, R.; Weihermüller, L.; Herbst, M. Future climate change will accelerate maize phenological development and increase yield in the Nemoral climate. *Sci. Total Environ.* **2021**, *784*, 147175. [CrossRef] [PubMed]

52. Herbst, M.; Pohlig, P.; Graf, A.; Weihermüller, L.; Schmidt, M.; Vanderborght, J.; Vereecken, H. Quantification of water stress induced within-field variability of carbon dioxide fluxes in a sugar beet stand. *Agric. For. Meteorol.* **2021**, *297*, 108242. [CrossRef]

53. Groh, J.; Diamantopoulos, E.; Duan, X.; Ewert, F.; Heinlein, F.; Herbst, M.; Holbak, M.; Kamali, B.; Kersebaum, K.C.; Kuhnert, M.; et al. Same soil, different climate: Crop model intercomparison on translocated lysimeters. *Vadose Zone J.* **2022**, *21*, e20202. [CrossRef]

54. Wallach, D.; Palosuo, T.; Thorburn, P.; Mielenz, H.; Buis, S.; Hochman, Z.; Gourdain, E.; Andrianasolo, F.; Dumont, B.; Ferrise, R.; et al. Proposal and extensive test of a calibration protocol for crop phenology models. *Agron. Sustain. Dev.* **2023**, *43*, 46. [CrossRef]

55. Brogi, C.; Huisman, J.A.; Weihermüller, L.; Klosterhalfen, A.; Montzka, C.; Reichenau, T.G.; Verrecken, H. Simulation of spatial variability in crop leaf area index and yield using agroecosystem modeling and geophysics-based quantitative soil information. *Vadose Zone J.* **2020**, *19*, e20009. [CrossRef]

56. Duan, Q.; Sorooshian, S.; Gupta, V. Effective and efficient global optimization for conceptual rainfall-runoff models. *Water Resour. Res.* **1992**, *28*, 1015–1031. [CrossRef]

57. Duan, Q.Y.; Gupta, V.K.; Sorooshian, S. Shuffled complex evolution approach for effective and efficient global minimization. *J. Optim. Theory Appl.* **1993**, *76*, 501–502. [CrossRef]

58. Wallach, D.; Palosuo, T.; Thorburn, P.; Hochman, Z.; Andrianasolo, F.; Asseng, S.; Basso, B.; Buis, S.; Crout, N.; Dumont, B.; et al. Multi-model evaluation of phenology prediction for wheat in Australia. *Agric. For. Meteorol.* **2021**, *298–299*, 108289. [CrossRef]

59. Wallach, D. Evaluating Crop Models. In *Working with Dynamic Crop Models*; Wallach, D., Makowski, D., Jones, J.W., Eds.; Elsevier: Amsterdam, The Netherlands, 2006; pp. 11–55.

60. Willmott, C.J. On the Validation of Models. *Phys. Geogr.* **1981**, *2*, 184–194. [CrossRef]

61. Coleman, K.; Prout, J.M.; Milne, A.E. RothC—A Model for the Turnover of Carbon in Soil. Model Description (February 2024). Rothamsted Research, Harpenden, 2024, p. 24. Available online: https://www.rothamsted.ac.uk/sites/default/files/Documents/RothC_description.pdf (accessed on 12 February 2025).

62. Vereecken, H.; Vanclooster, M.; Swerts, M.; Diels, J. Simulating water and nitrogen behaviour in soil cropped with winter wheat. *Fert. Res.* **1991**, *27*, 233–243. [CrossRef]
63. Aulakh, M.S.; Doran, J.W.; Mosier, A.R. Soil denitrification significance, measurements and effects of management. *Adv. Soil Sci.* **1992**, *18*, 2–42.

64. McIsaac, D.M.; Watts, D. *Users Guide to NITWAT—a Nitrogen and Water Management Model*; Agricultural Engineering Department, University of Nebraska: Lincoln, NE, USA, 1985.

65. Huwe, B.; Van der Ploeg, R. *Modelle zur Simulation des Stickstoffhaushaltes von Standorten mit Unterschiedlicher landwirtschaftlicher Nutzung*; Eigenverlag des Instituts für Wasserbau der Universität Stuttgart: Stuttgart, Germany, 1988; Volume 69, p. 213.

66. Feddes, R.A.; Kowalik, P.J.; Zaradny, H. *Simulation of Field Water Use and Crop Yield*; Simulation Monographs: Wageningen, The Netherlands, 1978; p. 188.

67. Vanclooster, M.; Viaene, P.; Diels, J.; Christiaens, K. *WAVE: A Mathematical Model for Simulating Water and Agrochemicals in the Soil and Vadose Environment: Reference and User's Manual (Release 2.0)*; Institute for Land and Water Management, Katholieke Universiteit: Leuven, Belgium, 1995; p. 154.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.