Nitrate vulnerability of karst aquifers and associated groundwater-dependent ecosystems in the Baltic region

Groundwater pollution by agrochemicals such as nitrogen fertilizers can cause complex biogeochemical transformations to take place in groundwater-dependent ecosystems. To explore the interaction between nitrogen load and groundwater-dependent, spring-fed ecosystems, a study was conducted in Latvia in an area of suspected high nitrate (NO3−) vulnerability due to its geological settings. A map of NO3− vulnerability along the margins of the carbonate aquifer in Latvia is presented. The map is based on a conceptual model that was developed during an extensive case study involving hydrological, hydrochemical, and habitat investigation of springs discharging from a karst aquifer and spring-fed ecosystems. Areas that should be prime targets for restricting fertilizer application are highlighted on the map. Although the case study revealed increased nitrogen pollution (up to 51 mg L−1, standard deviation of 9 mg L−1, in the springs discharging from the karst aquifer), no clear evidence of adverse effects due to NO3− pollution on the groundwater-dependent ecosystems using biotic indicators was found, highlighting the resilience of spring-fed ecosystems against high nitrogen inputs. In the case study, downstream groundwater-dependent ecosystems retained 70% of the reactive nitrogen during the vegetation season, but only a small proportion during the cold season. Thus, NO3− pollution can be partly mitigated by restoring wetlands along valley slopes where natural groundwater discharge takes place. The conceptual model developed for groundwater NO3− vulnerability is applicable to other areas in the Baltic region and other places with similar climatic and geological conditions.


Introduction
Groundwater pollution by agrochemicals such as nitrogen (N) fertilizers (Böhlke 2002) is a global problem that is causing the deterioration of the quality of water resources and is degrading groundwater-dependent ecosystems and other wetlands (Wulff et al. 2007;Gruber and Galloway 2008;Rockström et al. 2009;Giakoumakis et al. 2013). This is currently a problem in Latvia, where agricultural production has undergone significant intensification during the last 2 decades. The use of N fertilizers in Latvia has doubled from 2005 to 2019: from 40.9 to 80.7 thousand tonnes of N equivalent (Central Statistical Bureau of Latvia 2020). Inevitably, the increased use of N fertilizers has high potential to adversely affect groundwater quality and the health of groundwater-dependent ecosystems and other wetlands.
The N cycle is characterized by a complex biogeochemical transformation in the soil and groundwater. Generally, most N fertilizer leached into the groundwater is oxidized to nitrate (NO 3 − ) (Appelo and Postma 2005;Rivett et al. 2008). The main pathway of N attenuation in groundwater is bacterial reduction-denitrification (Rivett et al. 2008) to nitrite (NO 2 − ), then nitric oxide (NO), then nitrous oxide (N 2 O) and, finally, molecular N 2 . This process can take place in the absence of molecular oxygen, wherein electron donors are organic carbon or sulfide minerals (Böhlke 2002;Hansen et al. 2014). In favourable conditions, dissimilatory NO 3 − reduction to ammonia (NH 4 + ) (DNRA) also can take place (Necpalova et al. 2012). Excess nutrient inputs into natural ecosystems can lead to eutrophication and loss of biodiversity, but it can be triggered by excess SO 4 −2 inputs as well (Smolders et al. 2006;Geurts et al. 2009;Cirkel et al. 2014).
In the Baltic region, reducing conditions prevail in the aquifers (Højberg et al. 2017), favouring denitrification. However, the vast tile drainage network in the region bypasses the reducing groundwater system, and delivers the NO 3 − leaching from the soil into the surface waters (Højberg et al. 2017). In the past, observed NO 3 − concentration in groundwater were low (Levins and Gosk 2008) and sometimes elevated concentration of NH 4 + was observed (Klimas and Paukštys 1993). More recently elevated NO 3 − concentration (up to 50 mg L −1 ) was associated with shallow, often fissured aquifers associated with higher trace element concentration such as Cd, Mn, Ni, Pb, U and Zn (Retike et al., 2016a, b). Elsewhere, a range of contaminants have been found to be associated with NO 3 − pollution, delivered either directly by fertilizer input (As, Mn and Ni) or enhanced leaching from rocks and sediments (As, U, V, Cr) (Stamatis et al. 2011;Gamvroula et al. 2013).
Karst aquifers are particularly vulnerable to surface pollutants. They are characterized by triple porosity (matrix, fractures and conduits) that allows both laminar and turbulent flow (White 2002;Bakalowicz 2005;Worthington 2009). The duality of flow and the extreme heterogeneity and anisotropy of the permeability complicates the interpretation and sustainable management of karst aquifers (Ford and Williams 2007;Stevanovic 2015). The rapid transport and low retention capacity of pollutants in the karst aquifers renders them vulnerable yet important water resources (Bakalowicz 2005;Stevanović 2018) that can have seasonally variable water quality (Tsakiris et al. 2009;Alexakis and Tsakiris 2010), often supporting groundwater-dependent ecosystems (GDEs) at their discharge zones.
The European Union (EU) Nitrates Directive (91/676/ EEC) aims to reduce water pollution caused by N fertilizers used in agriculture by designating specific nitratevulnerable zones, where limitations on the application of N fertilizers and codes of good agricultural practice are implemented. The central part of Latvia was designated as a nitrate-vulnerable zone along administrative borders in 2001 (Cabinet of Ministers 2001), and it covers 12.8% of the area of Latvia. There is now political pressure to designate the whole country as a nitrate-vulnerable zone as previous studies of groundwater vulnerability have indicated that significant agricultural pressure exceeds the current boundaries of nitrate-vulnerable zone (Retike et al. 2016a). It has been shown that it is economically feasible to develop detailed denitrification maps and limit N fertilizer application only at sites that lack the natural capacity to cause NO 3 − reduction (Jacobsen and Hansen 2016). However, denitrification itself may cause undesirable consequences, such as mobilization of toxic heavy metals from oxidation of pyrite in sediments and release of nitrous oxide, a potent greenhouse gas (Palmer and Horn 2015). As the restrictions in the nitrate-vulnerable zone strongly interfere with agricultural practices, data-driven justifications are required to balance the needs for environmental protection and economic activity.
The aim of this paper is to present a map of groundwater nitrate vulnerability along the margins of a carbonate aquifer in Latvia. The map is based on a conceptual model that was developed during an extensive case study that involved hydrological, hydrochemical and habitat investigation of springs discharging from a karst aquifer and connected groundwater-dependent ecosystems. This conceptual model for groundwater nitrate vulnerability is applicable to the Baltic region and places with similar climatic and geological conditions. Furthermore, it was found that at such locations, fertilizers leaching from farmlands are partly retained by connected groundwaterdependent ecosystems. However, clear evidence of the adverse effects on ecosystems from such pressures was not found in this case study.

The Baltic region
The study site is located within Baltic Artesian Basin on the eastern shores of the Baltic Sea (Virbulis et al. 2013). The regional geology is characterized by a sub horizontal layering of Paleozoic and Mesozoic sedimentary sequences composed of both clastic and carbonate rocks. It is a few hundred metres thick in Northern Estonia and reaches several kilometres in thickness in Lithuania (Lukševičs et al. 2012). The land surface mostly was shaped by the Fennoscandian Ice Sheets and meltwater basins that formed during its retreat. The Quaternary cover ranges from a few metres in the plains to several hundred metres in some of the uplands, particularly in the southern and eastern part of the region. The terrain of the region comprises lowlands and hilly uplands that reach elevations of up to 318 m above sea level (a.s.l.). According to the Köppen-Geiger classification, the Baltic region has a warm summer, humid continental (Dfb, hemiboreal) climate (Klocking et al. 2006), but due to climate change starting from coastal areas of the Baltic Sea, it is shifting towards a temperate oceanic climate (Cfb).

The case study at Kazu Leja
Kazu Leja (also referred to as Kazu Grava) is a valley-like depression that is situated in Gauja National Park, a Natura 2000 site, in Latvia (Fig. 1,57° 19′ 58 N,25° 20′ 28° E). It is designated as a nature reserve zone and as a geological and geomorphological nature monument. The average annual average air temperature in the Priekuļi meteorological 2 km from the site is + 6.5 °C; January is on average the coldest month (-3.8 °C) and July the warmest month (+ 17.5 °C). Precipitation (711 mm year −1 ) is rather evenly spread across the year. The plains surrounding Kazu Leja have historically been intensively farmed arable land. Freshwater tufa and peat mining during the twentieth century left behind a heavily degraded landscape in the Kazu Leja itself. The contemporary habitats-fens, meadows and forestsare a result of ecosystem recovery. Numerous springs discharge along the slopes of the Kazu Leja. A particularly large group of springs is located in a side ravine (Lībāni-Jaunzemji ravine) near the western end of Kazu Leja (Fig. 1), feeding the Seven Springs Brook. A tile drainage ditch releases drainage water into the Seven Springs Brook as well. The spring water often infiltrates the scree deposits on the valley slopes and re-emerges at the surface in the valley floor as discrete springs or diffuse seeps. At the base of Kazu Leja, the Triečupīte River connects a string of fens, ponds left by peat extraction, and beaver-created channels, and merges with Seven Springs Brook before exiting the valley.
The Kazu Leja has a distinct U-shaped cross-section. It is 3.6 km long, 0.3-0.8 km wide and oriented E-W, connecting the valleys of the Vaive and Gauja rivers. The elevation of the surrounding undulating plain is 105-127 m a.s.l., and the valley floor is 65-72 m a.s.l. (Fig. 1). The plain surrounding Kazu Leja is covered with up to 10 m thick Quaternary till (diamicton, loam;Juškevičs 2000). The valley topography reflects a system of paleo-incisions in Paleozoic sedimentary rocks that was likely formed in proglacial or subglacial conditions during the last glaciations (Juškevičs 2000;Krievāns 2015) and partly filled with the layering of till and sand-gravel deposits (Bendrupe and Arharova 1981) (Fig. 2). The bedrock surface of the surrounding plain varies between 100 and 115 m a.s.l., but the base of the incised valley was determined by a borehole to be at about -30 m a.s.l. During the Holocene alongside the slopes and within the valley, peat and alluvial sediments and freshwater tufa accumulated (Āboltiņš 1998).
The bedrock is composed of an 8-13 m thick sequence of dolomites of the upper Devonian Pļaviņas Formation (D 3 pl) that overlie sandstones, siltstones, and clays of the upper Devonian Amata (D 2 am) and middle Devonian Gauja (D 2 gj), Burtnieks (D 2 br) and Arukila (D 2 ar) Formations. A clay-rich interval between the D 3 pl and D 3 am formations forms a regionally important aquitard (Virbulis et al. 2013). The upper part of the D 2 gj to D 2 ar formations are clay-rich intervals (Lukševičs et al. 2012), often considered to be locally porous aquitards, while the lower parts were considered as aquifers. Piezometric heads of these aquifers vary between 50 and 70 m a.s.l. (Fig. 2) (Spalviņš 2016). A 3D geological model of the Kazu Leja, developed as part of this study, is available at https:// www. gprm. lu. lv/ kazug rava-3d/.
Most springs in Kazu Leja emerge from the D 3 pl dolomite karst aquifer. Considerable amounts of freshwater tufa (Rozenšteins and Lancmanis 1924a;Melnalksnis et al. 1955) have precipitated from these springs during the Holocene (Timšs 1944). On warm summer mornings, occasional whitening events in the ponds collecting the spring water have been observed by locals, indicating ongoing carbonate precipitation. There are historical reports on sinkholes in the plains surrounding Kazu Leja (Rozenšteins and Lancmanis 1924b) indicating the presence of multiple karst features within the D 3 pl dolomites.

Methods
The research included extensive field campaigns at a pilot study site in Kazu Leja, examining habitats and hydrochemical and hydrological features as well as regional geological modelling. Data processing and computing were performed in the R environment (R Core Team 2020), particularly using the Tidyverse set of packages (Wickham et al. 2019).

Habitat mapping
Mapping of habitats was done to identify the groundwaterdependent ecosystems and to distinguish their types in July and August of 2019. The habitat types were determined according to the vegetation uniformity and abiotic conditions, and their relevance to the Latvian national interpretation (Auniņš et al. 2013) of the habitats listed in Annex I of EU Habitats Directive 92/43/EEC (Annex I habitat types) in 125 sample plots that were 50 × 50 cm large. In interpreting the ecology of plant communities (dominant and most frequent species), the Ellenberg indicator values from the Pladias database (Chytrý et al. 2021; https:// pladi as. cz/) were used.
The full methodology description is available at zenodo. org repository . In short, the major ions were analysed by liquid chromatography with SHIMADZU® RID-10A or CDD-10A conductivity detectors with Shodex IC YS-50 (cations) or Shodex IC SI-50 4E (anions) columns. Ion concentration was determined using multi-point calibration curves of the analytes prepared from standard solutions in the range of 0-100 mg/L (except 0-200 mg/l for Ca +2 ), with correlation coefficients of calibration curves r˃0.99. The DOC and N tot were analysed using Analytik Jena multi N/C 3100 TOC/TN b analyzer. P tot was analysed using standard methods (Murphy and Riley 1962;USEPA 1978;SFS 3025 1986;SFS 3026 1986) with Shimadzu UV-1800 spectrophotometer. Multi-point calibration curves of TC, TIC, N tot and P tot prepared from standard solutions in the range of 0-1000 mg/L (except 0-100 mg P L −1 for P tot ), with correlation coefficients of calibration curves r˃0.99. Trace elements in water samples were measured by Thermo Scientific Inc. ICP-OES spectrometer iCAP7000. Hard Drinking Water UK-Metals, ERM-CA011b, 0161 standard was used for measurement uncertainty quantification. Isotope ratios of hydrogen and oxygen were measured with Picarro L2120-i Isotopic Water Analyzer, a cavity ring-down laser spectroscopy method (Brand et al. 2009) following the IAEA protocols (Aggarwal et al. 2007) and expressed in δ-notation relative to the Vienna Standard Mean Ocean Water (VSMOW; Craig, 1961). Reproducibility less than ± 0.1‰ and ± 1‰ for δ 18 O and δ 2 H, respectively. The laboratory standards derived from IAEA standards (VSMOW2, GRESP) were used. The laboratory successfully participated in the IAEA isotope laboratory proficiency tests WICO in 2016 (Wassenaar et al. 2018) and 2020.
The water discharge was measured at locations Q1_railway, Q2_pond, Drainage_outlet and Spring_4) with salt dilution slug test introducing between 1 and 10 L of concentrated NaCl solution into the stream, measuring the EC at a downstream location with Diver®CDT logger. A calibration curve for each occasion was constructed to convert EC to salt concentration.

Reconstructing the subsurface catchment of the Kazu Leja
The subsurface catchments of the springs discharging into the Kazu Leja were estimated using a conceptual approach. As at the base of the D 3 pl aquifer, there is a subhorizontal aquitard. It was assumed that the groundwater in the D 3 pl aquifer would flow to the nearest location where this aquitard was outcropping at the sub-Quaternary surface. The distribution of the aquitard was modelled according to the methodology by Popovs et al. (2015). The location of the hypothetical groundwater water divide was estimated by constructing a Voronoi diagram from the nodal points of the linearly interpolated intersection between the surface of the aquitard and the land surface. A Voronoi diagram is a partition of the plane into regions closest to each of a given set of points. The Voronoi polygons attributed to the node-points within the surface catchment of Kazu Leja were considered to contribute their associated subsurface runoff to the springs of Kazu Leja.

Groundwater level in the D 3 pl aquifer
The groundwater table in the D 3 pl aquifer was estimated indirectly from the water table fluctuations in the inundated Lauciņi dolomite quarry near Kazu Leja. The water table elevation was estimated by deciphering the water line position of the quarry pond from aerial photographs at ten locations and reading the respective elevation on the digital elevation model constructed from the lidar data (LĢIA 2021).

Identification of sites with similar geological setting
Sites such as Kazu Leja, where NO 3 − contamination could be found in the D 3 pl aquifer in Latvia, were identified in a 0.25 km grid resolution based on the geological model. The modelled D 3 am surface (Popovs et al. 2015) that approximately coincides with an aquitard surface at the base of the D 3 pl aquifer was cropped by the sub-Quaternary surface elevation model. As the next step groundwater discharge sites from the D 3 pl aquifer identified were within a 0.75 by 0.75 km window, the lowest point of the land surface (LĢIA 2018) was below the D 3 am surface. Furthermore, the risk of NO 3 − contamination was assessed as a fraction of cultivated land from the CORINE land cover (European Environment Agency 2018) with a 4.75 km moving window accounting only for grid cells underlain by the D 3 am surface.

Habitats in Kazu Leja
During field campaigns in 2019, three terrestrial wetland and one aquatic habitat types were identified in Kazu Leja (Fig. 3); among them only one Annex I habitat type-Petrifying springs with tufa formation (7220*)-were found in Kazu Leja and had a total area of 2.3 ha. The mossdominated petrifying spring seepages and brooks with characteristic tufa-forming spring vegetation occurred in several locations covering patches of different size, forms few square metres to a relatively large complex with mixed slope forest with spring flushes, with an area of nearly 2 ha (Fig. 3). These contained a Cratoneurion commutati community with Palustriella commutata and rich in other bryophyte species. The most frequently recorded vascular plant species in these wetlands were Cardamine amara, Crepis paludosa, Myosotis palustris, Poa palustris, Carex remota). The fen on the valley bottom was not considered a protected habitat type according to the national habitat interpretation, as it does not host the typical plant assemblages, i.e. indicator species and their typical species associations of Annex I habitat types. However, in fact it is a spring-fed fen ecosystem that is species-poor, dominated by tall sedges (mainly Carex rostrata), Phragmites australis and Equisetum fluviatile vegetation. Besides the abovementioned plants, frequent and relatively abundant species in the fen included Epilobium palustre, Lysimachia thyrsiflora, L. vulgaris, Lycopus europaeus, Comarum palustre, with some patches of Carex paniculata and C. acutiformis. The fen spontaneously recovered after the cessation of peat cutting in the 1960s, which may be the main reason for lack of typical spring mire species associations. A part of the peatland on the eastern end of the valley has been afforested (Fig. 3). All wetland types in the valley (except at the eastern end of the area), both near-natural and those spontaneously recovering after severe damage in the past, depend on groundwater supply and function as wetlands, and thus were considered GDEs.

Subsurface catchment of Kazu Leja
The estimated area of the subsurface catchment of the D 3 pl aquifer that drains to Kazu Leja was 3.0 km 2 (Fig. 4), of which 1.8 km 2 was directed towards the Seven Springs Brook. The elevation of the reconstructed water table in the flooded Lauciņi dolomite quarry between 1997 and 2020 was estimated to be between 104.9 and 106.2 m a.s.l. Given that the elevation of many large springs in Kazu Leja is about 100 m a.s.l., and they are 1.8 km from the quarry, the linear groundwater gradient was between 0.0027 and 0.0034 m m −1 .

Discharge
The measured discharge of the Seven Springs Brook from the main group of springs was between 4.3 L s −1 in the summer (July) 2019 and 68.7 L s −1 in the early spring (March) 2020. Such a large variability reflects the flashy nature of the karst aquifer. During the field visits, the discharge from the agricultural drainage network was insignificant (up to 1.3 L s −1 ) or absent. However, locals reported occasional high flows in the drainage ditch associated with intense precipitation events, but the total runoff contribution from the field drainage seems small compared to the spring discharge. The average measured discharge of the Seven Springs Brook was 32.3 l s −1 , and thus the groundwater runoff would be about 566 mm year −1 .
The discharge of the Triečupīte River at its outflow from Kazu Leja was from 11.0 to 107.6 l s −1 . Of this, the Seven Springs Brook contributed more than 60% during peak flow (winter) and less than 40% during low flow (summer). The large springs appear to have a well-developed subsurface drainage network and therefore are more responsive to precipitation input (quick flow), while the rest of the valley delivers a greater share of the base flow.

Water chemical composition
This paper provides analytical results of 56 water samples collected from 11 locations in 2019 and 2020 in the Kazu Leja (Table 1). As summarized below, the dominant water type was Ca-Mg-HCO 3 and many of the samples had elevated concentration of biogenic elements (N and P).

Springs
The spring water had higher pH (between 7.1 and 8.4), EC (generally between 600 and 770 μS cm −1 ) and alkalinity (260 and 460 mg L −1 CaCO 3 ) values compared to samples that were collected from other sites. A slight downstream increase of pH was noted (Table 1), suggesting that CO 2 degassing was taking place. The molar ratios of Mg +2 to Ca +2 ions in Spring_2 and Spring_4 were 0.93 and 0.87, respectively. However, in other springs, the excess of Ca +2 over Mg +2 was slightly higher. Such observations are consistent with the almost stoichiometric weathering of dolomite in the D 3 pl aquifer.
The spring waters had a high level of dissolved O 2 (between 8.6 and 14.3 mg L −1 at temperatures between 7.0 and 8.6 °C), corresponding to low levels of dissolved organic carbon (up to 3.2 mg L −1 ) except for Spring_3 (up to 10.3 mg L −1 ). Spring_2 and Spring_4 was found to differ from the rest of the springs by having lower concentration of Na + , K + , Cl − and SO 4 −2 . The higher concentration of these ions in the main group of springs are likely due to the anthropogenic input of salts from fertilizers or road de-icing in the surrounding areas.
The stable isotope ratios fluctuated around the mean (δ 18 O − 10.5 to − 10.1‰ and δ 2 H − 74.9 to − 73.0‰; d-excess 7.5-9.5‰) aligned with local meteoric water line or was slightly below it (Fig. 5). It was suggested that intense precipitation leading to water impoundment on the soil surface can lead to groundwater isotope enrichment due to evaporation in the study region (Kalvāns et al. 2020). No distinct seasonal pattern of isotope ratios in spring water was noticed in contrast to the highly variable discharge. The median NO 3 − concentration in spring water was from 5 mg L −1 in Spring_2 to 33 mg L −1 (up to 51 mg L −1 ) in Q2_pond. The median concentration of total phosphorus (P) ranged between 0.02 mg L −1 in Spring_1 and 0.07 mg L −1 in Spring_4.

Agricultural field drainage
The pH of the drainage water varied between 7.5 and 8.1, the EC from 254 to 568 μS cm -1 and alkalinity between 120 and 240 mg-CaCO 3 L −1 . The median molar ratio of Mg +2 to Ca +2 in the drainage water was 0.41. This suggests that, unlike in the springs, the stoichiometric dissolution of the dolomite was not the main source of these ions. The relatively low EC levels are compatible with the fast runoff of the precipitation water.
The δ 18 O value of the drainage water ranged from − 10.4 to − 8.67‰, δ 2 H from − 72.8 to − 63.6‰, while d-excess from 5.7 to 11.3‰. The lowest d-excess value, indicative of surface evaporation, was observed in September 2019.
The concentration of NO 3 − in the drainage water was from 1 to 21 mg L −1 , while total P concentration ranged from 0.2 to 0.3 mg L −1 . The lowest concentration of NO 3 − in drainage water was observed in September 2019. The concentration of K + in the drainage water was very high (up to 9 mg L −1 ), and greatly exceeded the concentration measured at all other observation points (Table 1). This was probably due to K fertilizer application as it is one of the three major plant nutrients which rate is similar to P-fertilizers when expressed as pure oxides (Central Statistical Bureau of Latvia 2020).

Groundwater
The dissolved O 2 concentration in groundwater compared to other water types was low (median of 0.62 mg L −1 and 5.25 mg L −1 for Well_1 and Well_3, respectively); however, the O 2 measurement in Well_3 likely suffered from contamination during sampling due to water table dropping below the top of the screen interval. Generally, the sampled groundwater reflected the same geochemical traits as the spring water. The groundwater likely was derived from the re-infiltration of the spring water. However, the molar excess of Ca +2 over Mg +2 was slightly higher in the groundwater compared to the springs, possibly indicating the effects of the secondary dissolution of the freshwater tufa that is omnipresent in the topsoil of the Kazu Leja valley floor. However, concentration of Na + , K + , Cl − and SO 4 −2 observed in Well_1 were lower than in Well_3, but the opposite trend was observed in springs proximal to these wells (Spring_1 and Spring_2, respectively).

Table 1
Summary statistics of the water composition, see Fig. 1 for locations of the sites, full data set and methodology available at zenodo.org repository (Kalvāns et   Interestingly, a decreasing trend of both stable isotope ratios was observed with the highest values observed in spring 2019, and the lowest values in spring 2020. This might reflect the impact of the extraordinary drought in 2018. The NO 3 − concentration in the groundwater were moderate compared to the spring water, and the median was 4.1 mg L −1 and 6.9 mg L −1 for Well_1 and Well_3, respectively (Table 1). Well_3 had noticeably elevated concentration of NO 2 − (median 0.25 mg L −1 ), indicating ongoing NO 3 − denitrification (Rivett et al. 2008). Overall, the total P concentration in groundwater were similar to the upper range of the observed values in springs (median values of 0.05 mg L −1 and 0.09 mg L −1 for Well_1 and Well_3, respectively), except for one outlier from Well_3 (0.88 mg L −1 ).

Fen
Fen water was sampled in Well_2, where very high alkalinity (up to 680 mg-CaCO 3 L −1 CaCO 3 ) was associated with high Ca +2 concentration (up to 201 mg L −1 ). Secondary dissolution of tufa by acids produced by peat decomposition seems to be the main reason for the high alkalinity. In addition, the concentration of SO 4 −2 in these samples were lower than those of all other observation points (< 0.5 mg L −1 ), and the NO 3 − concentration was at most 0.67 mg L −1 . However, concentration of NH 4 + (up to 10.0 mg L −1 ) and total P (up to 0.34 mg-P L −1 ) was high. Apparently, the NO 3 − and SO 4 −2 were consumed by anaerobic peat decomposition, releasing some phosphorus; in addition, dissimilatory NO 3 − reduction to NH 4 − (Necpalova et al. 2012) also probably took place. Elevated concentration of Fe and Mn, which are active in reducing environments devoid of H 2 S, were noticed as well.
The fen water δ 18 O value fluctuated around − 9.4‰ and δ 2 H around − 69.6‰, and the d-excess value ranged from 5.3 to 6.6‰; seasonal variations were small, but the lowest d-excess was observed in July 2019. All samples had d-excess values indicative of evaporation (Fig. 5), suggesting that the fen water residence time is greater than 1 year.

Triečupīte river
The most notable characteristic of the river water (River in Fig. 1) was its higher DOC concentration (from 2.4 to 32.3 mg-C L −1 ) than at other locations ( Table 1). The highest DOC concentration at this site was observed in July 2019 indicating high biological activity and perhaps eutrophic conditions and more intensive leaching of organic substances from the peat substrate. The river water had relatively high alkalinity (ranging from 280 to 460 mg-CaCO 3 L −1 ). High alkalinity and DOC could be the result of internal eutrophication facilitated by the high buffering capacity of omnipresent freshwater lime (Smolders et al. 2006).
The molar ratio of Mg +2 to Ca +2 varied with a median of 0.80, apparently inherited from upstream spring water. The low concentration of Na + , Cl − and SO 4 −2 was noteworthy, even lower than those observed in Spring_2 and Spring_4. This might be due to dilution with a precipitation input to the catchment part other than from the carbonate aquifer.
The δ 18 O value was from − 10.8 to − 8.0‰, δ 2 H from to − 76.5 to − 62.7‰ and d-excess from 1.4 to 10.5‰ in the river water. Very clear seasonality was observed, with most enriched ratios observed in summer 2019, plotting below the meteoric water line in the dual isotope plot (Fig. 5) and indicating isotope enrichment due to evaporation. NO 3 − concentration in the river water were negligible in the warm season (May, July, September 2019), up to 1.1 mg L −1 , and rather high in colder months (November 2019 and March 2020), 12.0 and 19.8 mg L −1 . Such a seasonal pattern points to the consumption of NO 3 − by biological activity. Total P concentration did not show such a clear seasonal behaviour as they ranged from 0.05 to 0.30 mg-P L −1 , with the highest concentration being observed in July and November 2019. On average, total P concentration in the river were higher compared to the springs, but NO 3 − concentration was smaller than in most of the springs.

Mixed stream
Water sampled near the confluence of the Triečupīte River and Seven Springs Brook (sampling point Q1_railway) had a chemical composition that was between the spring water, particularly Spring_3, and the river water. Its pH was rather stable around 7.9 while EC and alkalinity values were similar to Spring_1 (640 μS cm −1 and 300 mg-CaCO 3 L −1 ). Likewise, the molar ratio of Mg +2 to Ca +2 of the mixed water was narrowly confined and like Spring_1, i.e. between 0.60 and 0.75. The stable isotope ratios for mixed water had a seasonal character inherited from the river water but with a smaller amplitude.
The total P concentration (0.05 mg L −1 on average) was higher while total N (average total N was 5.7 mg-N L −1 ) lower than in most spring water samples but lower than that of the river samples.

Sites with similar geological and land use conditions
A geological model was used to identify locations with similar settings as the Kazu Leja case study site. These were located along the Gauja River and other river valleys close to the distribution margin of the D 3 pl aquifer (Fig. 6). The relevant GIS data layer is available at zenodo.org (Kalvāns and Popovs 2021). The modelled discharge zones and actual spring locations in Kazu Leja showed reasonable agreement (Fig. 4) given the coarse resolution of the model (250 m).

Discussion
A map of groundwater NO 3 − vulnerability for the D 3 pl aquifer in Latvia was developed. It is based on a conceptual hydrogeological model for the Kazu Leja-a valley cut into plain terrain-case study. According to the conceptual model, most of the water reaching Kazu Leja discharges from springs that drain the D 3 pl aquifer. The large seasonal fluctuations of spring discharges from the D 3 pl aquifer, as well as historical evidence of the presence of sinkholes, indicate that this is a karst aquifer. The spring water had elevated NO 3 − concentration, but the downstream fen retained some of the excess nitrogen. The surface runoff into Kazu Leja seems to be insignificant compared to the natural groundwater discharge. The groundwater discharge is crucial for the recovery and sustenance of historically degraded petrifying spring and fen ecosystems on the valley slopes and floor. It is likely that the conceptual model of groundwater flow and NO 3 − pollution is applicable to other geologically similar locations in the Baltic region and around the world.

Geological conditions were favourable to NO 3 − conservation in the groundwater at the Kazu Leja site
The high NO 3 − concentration in the springs can be explained by geological conditions-a thin layer of Quaternary sediments covering the D 3 pl karst aquifer are favourable to NO 3 − conservation in groundwater due to the well-aerated state of the aquifer (Fig. 7). The vadose zone remains within the dolomite aquifer for at least several kilometres from the spring discharge sites, as seen in the flooded Lauciņi dolomite quarry; thus, the infiltration water remains oxygenated, preventing NO 3 − leaching from the soil. The Devonian sediments in the study region contain only refractory organic substances, which are not a suitable substrate for microbiological activity (Prols 2010) and do not allow for microbial NO 3 − reduction. Similar patterns have been observed elsewhere; for example, investigating a karst aquifer covered by thin Quaternary cover in north Estonia, (Koit et al. 2020) found that carbonate dissolution took place in open system conditions, implying that an aerated environment provided conditions favourable to NO 3 − conservation with the vadose zone. Indeed, elevated average NO 3 − concentration (up to 23.2 mg L −1 ) were found at this location (Koit et al. 2020). Thus, the conceptual model of NO 3 − vulnerability of karst aquifers appears to be applicable around the region.
Considering the catchment area and discharges of the springs, the average intensity of NO 3 − runoff from the D 3 pl aquifer was 50 kg N ha −1 year −1 (Table 2), close to the country-wide average use of N fertilizer in Latvia during 2015-2019 (64 kg-N ha −1 year −1 (Central Statistical Bureau of Latvia 2020). However, the calculated groundwater runoff (566 mm year −1 ), and thus NO 3 − leaching rate, seems to be overestimated compared to the long-term average runoff for the Gauja River basin (247 mm year −1 ; Klavins et al. 2002). It is likely that the recharge area for the springs was underestimated, explaining the overestimation of the runoff. As mineral soils and farmland predominate in the catchment area, Fig. 7 A conceptual model of groundwater recharge and discharge favouring NO 3 − conservation, where a thick aerated zone was formed in a karst aquifer near a cliff-like incision, and where NO 3 − retention takes place in a connected groundwater-dependent ecosystem the NO 3 − pollution must come from recently applied fertilizers, rather than legacy pollution or peat decomposition.

Seasonally favourable conditions for NO 3 − retention
The results indicate that a noticeable retention of reactive N takes place in the wetland at the base of the Kazu Leja during the warmer season of the year. The N concentration in the outflowing water-the Triečupīte River-was negligible in the warm season, but comparable to springs during the cold season. The average total nitrogen concentration in the Triečupīte River above its confluence with Seven Springs Brook in summer was three times lower than in winter (1.13 and 3.37 mg-N L −1 , respectively) and the summer concentration of NO 3 − was negligible. Thus, assuming an insignificant seasonal pattern of N input, it can be calculated that total N retention was about 70% during the warm season, but only 15% on average for the year.
The local topography-a flat valley floor and plentiful groundwater supply-favours the establishment of hygrophilous vegetation and the development of GDEs capable of NO 3 − retention. Indeed, artificial wetlands are used for wastewater treatment, but are known to perform poorly in cold climates, during winter (Varma et al. 2021) when microbial activity decreases and macrophyte growth ceases. Thus, it is tempting to consider the GDEs at such sites as natural water filtration facilities. However, the denitrification itself can have undesirable consequences, such as mobilization of toxic heavy metals from oxidation of pyrite in sediments and release of nitrous oxide-a potent greenhouse gas. Indeed, elevated concentration of heavy metals associated with agricultural pollution have been found in Latvia (Retiķe et al. 2016a, b).

Vegetation response to high nutrient input by water
The P concentration in the Triečupīte River and fen water significantly exceeded the background concentration in potomal small river in the Gauja River basin (0.045 mg-P L −1 ; LVĢMC 2015), but was rather moderate in the springs. In contrast, the total N background concentration (1.5 mg-N L −1 , LVĢMC 2015) was similar to background but significantly exceeded the springs. Thus, it can be concluded that the P source in the river was the wetlands itself rather than external input.
In Kazu Leja, the SO 4 −2 was absent in the fen water (Table 1, Well_2), while an elevated concentration of Fe and elevated concentration of P were observed. The internal eutrophication of wetlands can be triggered by excess input of SO 4 −2 that reduces to H 2 S, which can be toxic to higher plants (Geurts et al. 2009) and bind divalent iron into sulfide minerals. As a result, the bivalent iron is not available to form trivalent iron minerals that bind phosphate ions where favourable redox conditions prevail (Smolders et al. 2006). Thus, it can be inferred that at the sampling location, any H2S produced by SO 4 −2 reduction was precipitated as iron sulfides by an excess of bivalent Fe. Whether this pathway of eutrophication was important in Kazu Leja remains to be clarified by further investigation.
There was no evidence of unfavourable impact on the fen ecosystem due to increased nutrient input, i.e. no clear signs of fen eutrophication using plants as indicators. According to the Ellenberg indicator values, some of the most dominant and common species found in the Kazu Leja fen (Carex rostrata, Epilobium palustre, Lysimachia thyrsiflora) indicate poor to intermediate nutrient conditions, but very rarely occurring at nutrient-rich sites. Some typical species (Equisetum fluviatile, Carex paniculata, Lycopus europaeus) prefer moderately rich sites, rarely occurring in extremely rich or poor sites. However, only a few typical species in the Kazu Leja fen (Comarum palustre and Carex diandra) are found in nutrient-poor conditions. No species affiliated with narrow ecological niches, e.g. clearly related to solely nutrient-rich or nutrient-poor conditions, were present. Overall, the plant community is typical for a moderately rich fen. The land use in the spring catchment area has not changed for several decades, and thus the current plant species composition is the long-term result of interaction among various environmental drivers and biota. The results suggest that both the fen and petrifying spring plant communities are resilient regarding increased nitrogen load.

Similar geological conditions identified along periphery of the D 3 pl aquifer
A map of locations shows where groundwater contamination with NO 3 − from arable land was likely (Fig. 6) and ), although the presented study was limited to the D 3 pl aquifer in Latvia. However, there are other equally significant aquifer formations in the Baltic region. The famous Estonian Klint-an escarpment up to 56 m high and several hundred kilometres long (Raukas 1997) largely composed of Ordovician and Silurian carbonate rocks (Perens and Vallner 1997)-is a prime example where such conditions might prevail. The Estonian Pandivere upland, composed of carbonate rocks (Vallner 1996) and surrounded by extensive wetlands, is another example. In Latvia, there are a number of Upper Devonian carbonate formations exposed at the bedrock surface and cut into by numerous river valleys (Lukševičs et al. 2012). The presented conceptual model can be adapted for NO 3 − vulnerability mapping in these regions as well.
The karst aquifer that discharge water to the Kazu Leja has no capacity to reduce NO 3 − pollution; however, the connected wetland ecosystem does. Delineation of the sites where geological conditions are not favourable to N 628 Page 14 of 16 retention in groundwater can be useful both for estimating the groundwater contribution to NO 3 − retention (Højberg et al. 2017) and for designing economically sensible guidelines for limiting the fertilizer application (Jacobsen and Hansen 2016). On one hand, in such locations, natural or constructed wetlands may play a crucial role in nitrogen retention as a nature-based solution. On the other hand, increased nutrient concentration may be detrimental for fen species and communities in nutrient-poor conditions, causing loss of threatened or otherwise important species and turnover of species composition (Smolders et al. 2006;Cirkel et al. 2014), and can increase greenhouse gas emissions. In the case of Kazu Leja, no direct evidence that the elevated nutrient input had a profound negative impact on the ecosystem was found, but its status as a nature reserve zone of the Gauja national park indicates that priority should be given to restoration or self-recovery of the fen ecosystem. In similar sites without considerable biodiversity potential, their nutrient retention capacity could be considered as a part of an integrated water resources protection plan.
In this study, one scenario illustrating the interaction between land use and groundwater systems was highlighted. The management actions to be undertaken form an open question that deserves further investigation, including the involvement of stakeholders.

Conclusions
A map indicating locations where the D 3 pl karst aquifer in Latvia is sensitive to nitrate (NO 3 − ) pollution was presented. The map was based on a conceptual hydrogeological model of a case study in Kazu Leja, a groundwater-dependent ecosystem. During the case study, elevated NO 3 − concentration (up to 51 mg L −1 ) in spring water were found. Geological conditions-thin Quaternary cover overlying unconfined karst aquifer-rendered the groundwater sensitive to the NO 3 − pollution. The downstream wetland ecosystem retained about 70% of the reactive nitrogen during the warm season but an insignificant amount during the cold season. The elevated nutrient input had no profound negative impact on the ecosystem. Therefore, it is concluded that wetland conservation, restoration or construction of artificial wetlands could mitigate the downstream export of the NO 3 − contamination. However, farmland runoff can lead to undesirable processes in ecosystems, that were not detected in this study. The results can inform stakeholder discussions and planning of water resources management.