Streamflow–concentration relationships of surface water in the Choapa basin: historical analysis and projections under climate change

ABSTRACT The use of concentration–streamflow (C–Q) relationships is presented as a novel approach for studying climate change effects on water quality. Based on data from nine monitoring stations in the Choapa basin, north-central Chile, constituent behaviours were classified as constancy, enrichment or dilution. Constancy was shown for B, As, and Cu. Dilution and enrichment relationships were observed at some sites. Electrical conductivity and major ions showed dilution behaviour, with Na+, Mg2+, Cl−, SO4 −2 and K+ showing greater variability in response to streamflow changes. Fe, Al, Mn and Zn presented enrichment behaviour, with Al and Fe showing greater variability. Based on historical C-Q relationships and available projections of streamflow variations under climate change, the change in constituent concentrations from 2010 to 2040 likely will not exceed ±10% with respect to the historical average, and ±15% from 2040 to 2070. In particular, Fe and Mn require special attention in the future.


Introduction
Climate change potentially affects both water availability and water quality. Since 2010, Chile has experienced a period of extreme water scarcity, known as the "megadrought" (Núñez andVerbist 2018, Garreaud et al. 2019). This has resulted in a deficit of approximately 70% in average streamflow in regions such as Coquimbo and Valparaíso (north-central Chile), as well as a continual decrease in groundwater levels (Valois et al. 2020). A total of 82 communities have been considered "under water scarcity" between the Atacama and Aysén regions in the period 2008-2017, mainly concentrated in the Coquimbo, Valparaíso, Maule and Metropolitan regions, between 30 and 34°S (MMA 2018). After 2021, the entire territory of the Coquimbo Region has experienced reduced streamflow in all rivers due to precipitation deficits (MOP 2021).
The country includes at least seven of the nine "vulnerable landscapes and characteristics" traits defined by the United Nations Framework Convention on Climate Change (Muck 2012, MMA 2016, Barrera et al. 2020. Therefore, in addition to the critical current situation, future hydrological projections in Chile are rather worrisome. These projections include (a) a decrease in precipitation (by 2030) of 5-15% in a large part of the territory (between 26 and 49°S), which would intensify by 2050; (b) a marked increase in drought events, especially from the second half of the 21st century; and c) a significant decrease in average monthly flows, again in a large part of the territory (30-42°S) with an elevation in the 0° isotherm, and a retreat of glaciers (MMA 2017).
Besides water quantity-related forecasts, there are potential climate change effects on water quality, but these are less well understood and have received less attention in both the international scientific literature and in Chile (Byrne et al. 2020, Peña-Guerrero et al. 2020. This paper therefore contributes appropriate approaches to address effects on water quality, specifically for data-limited regions where complex (and datademanding) modelling efforts are not feasible (Nauditt et al. 2017a). Within this context, we explore the study of historical relationships between concentrations (C) of various water constituents and streamflow (Q) -in short, the C-Q relationships (Bonta 2005). These relationships "represent the integration of biogeochemical and physical processes within the hydrologic source areas that contribute to streamflow" (Rose et al. 2018(Rose et al. , p. 2830. Thus, and despite the inherent uncertainties about future climate conditions (e.g. Souvignet et al. 2010), proper knowledge of the C-Q relationships for different water constituents would allow us, at least in a preliminary stage, to extend available climate change forecasts for future hydrological conditions to water quality issues.
Although different types of C-Q relationships have been described, the one most commonly used is the power law relationship of the form C = aQ b (Flores et al. 2016, Li et al. 2021. This relationship allows classification of water quality conditions as chemostatic (constancy), i.e. concentration independent with respect to streamflow, or non chemostatic. The latter can be divided into two categories: (a) C increases as Q increases, referred to as enrichment, flushing or enhanced hydrological access (Salmon et al. 2001, Li 2019; and (b) C decreases as Q increases, referred to as dilution. Thompson et al. (2011) and Musolff et al. (2015) identified cases where the slope "b" is 0 (or not statistically different from 0), with high C variability relative to Q, which they denominated as chemodynamic (see also Pohle et al. 2021). Briefly, regarding specific behaviours for particular water constituents, Godsey et al. (2009) analysed the behaviour of different elements in 59 US watersheds, identifying that Ca 2+ , Mg 2+ , Na + , and Si generally showed chemostatic behaviour. Later, Hunsaker and Johnson (2017) analysed C-Q relationships in eight experimental watersheds in the King River area of southern Sierra Nevada, California (USA), from 2004 to 2011. These authors found that Ca 2+ , Mg 2+ , K + , Na + and Cl − presented negative slopes ("b" values), close to but statistically different from 0, differing from Godsey et al. (2009). Additionally, Rose et al. (2018)  , total phosphorus, and total suspended solids showed increases in C with increased Q. Likewise, Zhu et al. (2018), identified heterogeneous behaviours for trace elements (e.g. Cu, Fe, Co, Ni, As, among others), in response to variation of Q in the Min River, Tibetan Plateau area. For the sake of space, we refer to Musolff et al. (2015), Li (2019), and Pohle et al. (2021) for a more in-depth analysis of the nature and causes of constancy (chemostaticchemodynamic), enrichment or dilution patterns.
In Chile, concentration and streamflow have shown both inverse and direct relationships. In the Choapa River basin (southern part of the Coquimbo Region, north-central Chile), Parra et al. (2011) qualitatively described an inverse relationship of As concentrations with Q in the upper zone of the basin, while Fe showed an increase in C with increasing Q, and Cu showed less dependence. In the lower areas of the basin, they observed an inverse C-Q relationship for the major anion concentrations (HCO 3 − , SO 4 2− and Cl − ). In the Andean tributaries of the Elqui basin, in the northern part of the Coquimbo Region, Flores et al. (2016) found an inverse semi-logarithmic relationship between electrical conductivity (EC), Ca 2+ , Mg 2+ , and SO 4 2− streamflow, while only Fe showed a mostly direct relationship. Based on these results, a very preliminary (and qualitative) assessment of the possible effects of climate change on water composition was carried out in the upper part of this basin, identifying rather favourable forecasts for Fe, but the opposite for SO 4 2− , in that a future decrease in flows would result in lower and higher concentrations of those constituents, respectively. More recently, Peña-Guerrero et al. (2020), using historical data from the Maipo River basin (central Chile), identified that EC in irrigation water was the parameter most influenced by drought, increasing as flows decreased, even exceeding the Chilean standard for irrigation water use, NCh 1333 (INN 1987). They also pointed out that water quality depended on several factors, one of the main ones being anthropogenic intervention, thus explaining the behaviour of variables that were not directly related to streamflow.
Beyond the specific studies referred to in the previous paragraphs, and to the best of our knowledge, more complete and advanced analyses have not been developed in Chile regarding C-Q relationships, nor about the possible effect of climate change on water quality, at least in the central and central-northern areas of the country, i.e. between 29 and 35° S. This is a concern considering that more than 75% of the country's agricultural activity and some of its most populated cities are concentrated in that area (Demaria et al. 2013, Garreaud et al. 2017). In addition, several watersheds in this area, particularly in headwater zones, present mining activities and natural conditions, e.g. presence of hydrothermal alteration zones (Zegers et al. 2021, Montserrat et al. 2022) and rock glaciers (Schaffer et al. 2019). Thus, climate change effects on the hydro-meteorological regimes and their relation to water quality in this region are still poorly understood, despite the inherent importance of this issue. This is particularly relevant as mountainous headwater basins have been recognized in recent years to be of great importance worldwide, as they constitute true "water towers" by accumulating, regulating, and delivering water to downstream areas (Viviroli et al. 2007, Nauditt et al. 2017b. Also, high-elevation mountains worldwide have been identified as areas highly sensitive to climate change (Li et al. 2021), as there is strong evidence that the rate of warming increases with elevation (Zarroca et al. 2021).
The present work presents a novel approach for the use of C-Q relationships in climate change-related water quality studies in mountain river systems of semi-arid zones. A case in point is the Choapa River basin, which is the subject of the present work. The Choapa River basin is a mountainous river basin with a short distance from the headwater to the river outlet, which features agriculture and mining activities. As such, this basin can serve as a valuable natural laboratory for the study of the C-Q relationship for different water constituents and under different water-use activities. Thus, this paper has two main objectives: (1) to characterize C-Q relationships of constituents in the surface waters of the Choapa River basin; and (2) to infer, from these relationships, the effects of forecasted climate change on streamflow and the behaviour of the water constituents. Beyond the results for the Choapa basin, the methodology proposed can be of interest and use for similar basins elsewhere, specifically with activities such as agriculture and mining, and generally when subjected to the effects of climate change on hydrology and river water quality.

Study area
The Choapa River basin ( Fig. 1) is located in the Coquimbo Region in north-central Chile, between 31°10' and 32°15'S, in the narrowest part of the country (Paskoff 1993), covering an area of 8124 km 2 . The Choapa River originates in the Andes Mountains at the confluence of the Totoral, Leiva and del Valle rivers, about 140 km east from the Pacific Ocean. Its initial, main tributaries are the Cuncumén and Chalinga rivers, the first of which receives drainage from the Los Pelambres (MLP) mining area. In the centre part of the basin, the Choapa River joins the Camisas stream and then joins an important tributary of the basin, the Illapel River. In the lower reaches it meets the Canela and Millahue streams (DGA 2004).
Precipitation is concentrated in the winter months (May-August), reaching a total of 270 mm/year (long-term average). Maximum annual temperatures occur in the summer months (December-March) and minimum temperatures in winter (May-August). Streamflows show a characteristic response of a snowy regime, with maximum runoff observed in the months of November and December (Vicuña et al. 2020).
Rock glaciers in the Choapa River basin comprise a total area of 18.4 km 2 , in the headwaters upstream of Salamanca, while no debris-free glaciers are present Brenning 2010, Schaffer et al. 2019). Rock glaciers in the catchment are relatively small, and are somewhat degraded, due to both environmental conditions and the impact from mining in the area (Monnier and Kinnard 2013). Whilst there is a possibility that they contribute to their local catchment, the likelihood of a significant contribution to downstream flow, compared with precipitation, is constrained by their relatively small areal coverage (Schaffer et al. 2019).
The upper zone of the basin is dominated by volcanosedimentary rocks, the middle-upper zone by stratified dioritic to granodioritic granitoids that host Fe, Cu and Au mineralization, the central zone by marine and transitional sedimentary sequences of sandstones, conglomerates, and limestones, and the lower part of the basin by granodiorites, tonalites and diorites, as well as monzodiorites, gabbros and diorites. The geological formations on which the riverbeds are developed correspond to unconsolidated detrital fillings of alluvial origin, mainly gravels, sands and silts from the rivers or their terraces and floodplains (SERNAGEOMIN 2003, DGA 2004. Two main economic activities, agriculture and mining, are carried out in the basin, with agriculture using the greatest amount of water (DGA 2017). Based on data from (Laboratorio de Prospección, Monitoreo y Modelación de Recursos Agrícolas y Ambientales, PROMMRA, w/o year), in 2016-2017, a total of ca. 9100 ha were cultivated in Choapa, mainly with fruit trees and annual crops. Annual average water demand for agriculture is on the order of 9800 L/s (canal intake level). This high water demand is due to the use of lowefficiency irrigation methods in most (ca. 77%) of the cultivated area.
Mining in the headwaters of the basin includes the worldclass porphyry copper deposit of MLP (Parra et al. 2011). Also, significant mining activity on a smaller scale (both active and abandoned) is scattered throughout the area (Fig. 2).
Two main cities in the basin are Illapel and Salamanca, home to about 60 000 people. A number of smaller settlements are also distributed throughout the basin, which are served by more than 30 rural drinking water systems (SSR in the Spanish acronym). Two main water reservoirs are present in the basin: the Corrales (50 Mm 3 , since 2000) and El Bato (26 Mm 3 , since 2012).

Database
The historical series of daily averaged Q (streamflow) and C (concentrations for solutes, or levels for pH or EC) were obtained from the General Water Authority (Dirección General de Aguas, DGA), available from (Sistema Nacional de Información del Agua, SNIA) (https://snia.mop.gob.cl/ BNAConsultas/reportes). Table 1 presents characteristics of the stations (for location see Fig. 1) and the extent of the original data series. While streamflow data are generally available at a daily time scale, water sampling by the DGA (for chemical composition) was conducted quarterly over much of the extent of the record, and bimonthly since about 2005. The DGA water constituent data correspond to total concentrations, based on sampling and analytical procedures described in the "Standard Methods for the Examination of Water and Wastewater" (Peña-Guerrero et al. 2020). For details of the analytical methods followed by DGA, the reader is referred to Oyarzun et al. (2006) and Flores et al. (2016).
Given the characteristics and objectives of this study, Q and C data (for the different hydrochemical variables) were concurrent, i.e. generally for the same day. However, for 13 cases (of the entire dataset considered), when the DGA database had missing flow data for a specific day corresponding to concentration data, Q of the closest available day was used. This gap between the date of water quality sampling and available streamflow data was only one day in 10 of the 13 cases identified, and never exceeded five days (the latter occurring on only one occasion).
Once the initial database was established, values with "lower than detection limit" were eliminated from the analysis (i.e. C-Q relationships). Subsequently, an identification of outliers was performed. Two complementary approaches were considered: (1) box plots, which were created with the "graphics" package in R, version 4.0.2 (R Core Team 2014); and (2) Test Box B, included in the "univOutl" package of R (D'Orazio 2018). Outliers were individually evaluated to determine whether they could be considered valid records, based on the following criteria: (a) comparison with the behaviour of other parameters on the same date (e.g. major ion and/or EC); (b) comparison with the same parameter at another nearby monitoring station on similar dates; and (c) ionic balance of the sample. Unreliable outliers were discarded from the database as they most probably corresponded to an error in the sampling/analysis and/or the data transcription. A final criterion for including a parameter in the analysis was that a minimum number of 30 C-Q pairs were available. Table 2 shows the resulting availability of data for each station and parameter.
With respect to the parameters selected, there are different classifications in the literature. For example, Rose et al. (2018) classified them as geogenic or weathering products (Ca 2+ , Mg 2 + , Na + , Si), exogenous solutes (Cl − , NO 3 − , SO 4 2− ) and biologically associated parameters (NH 4 + , PO 4 3− , Dissolved Organic Carbon (DOC), Dissolved Organic Nitrogen (DON), K + ), whereas Zhu et al. (2018) considered a classification of major elements and trace elements. For the sake of clarity and simplicity, in the current work, we followed Musolff et al. (2015), so we considered a simple classification into field parameters, major ions, nutrients, and trace elements.

Descriptive analysis and C-Q relationships
Initial descriptive analyses considered time series scatter plots and box plots for each constituent and station. Then, simple regressions were determined using the R software ("stats" package), based on the historical values of both Q and C data for each station and for each parameter (in log-log plots), resulting in parameters "a" and "b" in the relation C = aQ b . The C-Q relationships were classified depending on the value of the regression slope "b." The literature does not suggest a single criterion (i.e. threshold value for "b") to use in classifying the behaviour of a particular parameter or water constituent. Godsey et al. (2009) simply stated that a slope of 0 represents a chemostatic behaviour. Li et al. (2017) defined chemostatic behaviour as occurring when −0.2 < b < +0.2, whereas Li (2019) and Ackerer et al. (2020) considered that this behaviour can be identified by a condition of −0.1 < b < +0.1. The latter was also assumed by Herndon et al. (2015), without considering as a requirement that the slope is statistically different from 0. More recently, Pohle et al. (2021) differentiated between enrichment (positive "b"), dilution (negative "b") and constancy ("b" ~ 0), determined though Mann-Kendall testing, and chemostatic/chemodynamic behaviours. Following Hunsaker and Johnson (2017), instead of the absolute value of the slope, we used a t-test to determine whether "b" was significantly different from 0. The following classification was followed in this work: (a) constancy behaviour was considered when "b" was not significantly different from 0, at a significance level of 5%; (b) when "b" was significantly different from 0, the C-Q relationship was classified as dilution if the "b" was negative, and as enrichment if "b" was positive.
We also calculated the CV C /CV Q index (i.e. the relationship between the coefficient of variation of C and Q) for each parameter and monitoring station, following Musolff et al. (2015), Rose et al. (2018) and Zhu et al. (2018). The advantage of considering CV C /CV Q , in complement to "b," is that the index is not affected by the initial assumption of a power law relationship (i.e. C = aQ b ). Thus, the combined use of the C-Q slope (b) and CV C /CV Q provides a more robust approach to characterize the behaviour of a given water constituent (Rose et al. 2018). Regarding this index, we followed Musolff et al. (2015) in that a CV C /CV Q value close to or greater than 1) represents a high degree of C variability with respect to Q, i.e. a chemodynamic behaviour, whereas CV C /CV Q ≪ 1 is representative of a chemostatic behaviour. The characterization of C-Q relationships allowed us to estimate possible effects of climate change on water quality, as described next.

Climate change effects
Regarding the second objective of this work, we evaluated the possible effects of hydrological changes, as a consequence of climate change, on water composition of the rivers in the study area.
We considered available studies that present predictions of precipitation and temperature, and from them, of future Q values, for the Choapa basin. These studies were: (1) CCG (2013), which considered 15 global circulation models and two GHG emission scenarios (A1B and B1); and (2) Vicuña et al. (2020), which was based on six climate models and the Representative Concentration Pathways (RCP) 8.5 highemission scenario. Complementarily, Montserrat et al. (2022) presented estimates of Q changes for the nearby upper Mapocho River basin, based on 10 global circulation models and the RCP 8.5 scenario, which was also considered. We are confident that these studies provided representative Q estimates because they reflect the hydrology of the basin where the focus on discharge is mostly determined by snow accumulation and melt processes in Andean Pre-Cordillera and Cordillera areas. From the mentioned works, we considered reductions of 14% and 20% in Q (spatially averaged annual flow rates) for the periods 2010-2040 and 2040-2070, respectively. Moreover, these projected (i.e. climate changerelated) Q values were evaluated to ensure that they were within the historical range of the observed data (C-Q relationships) before estimating future water constituent concentrations. In this way, the regression parameters were used for the purpose of interpolating (not extrapolating) the Table 2. Number of validated records of streamflow-concentration data pairs ("-" denotes parameters not considered at a given station due to insufficient data). Class  E1  E2  E3  E4  E5  E6  E7  E8  E9  Field parameter  pH  281  104  62  209  94  148  118  118  232  EC  281  104  62  209  94  148  118  118  232  Major ion  Na +  182  92  60  160  58  107  61  109  178  Ca 2+  183  94  61  162  58  108  61  110  179  Mg 2+  182  94  61  161  58  108  61  110

Parameter Station
estimated concentrations for a given future streamflow. Also, estimates of change in future concentrations were only made for those parameters and stations for which dilution or enrichment behaviours were identified (based on the C-Q relationships obtained from historical data). Using the C-Q relationships established from the historical database, we calculated the variation in C (ΔC, in %) due to projected variations in Q (ΔQ, in %), as a function of regression slope (b), with the equation ΔC = (1 -ΔQ) b -1 (Table 3). These projected variations were applied to the historical averages of the constituents selected. The rationale for using averages and not medians was to perform an analysis in a more unfavourable condition, as averages are normally greater than medians, especially for trace elements (e.g. Cu, Fe, As).
Finally, we compared these projected C values to the limits established by NCh 1333 for irrigation (INN 1987) and NCh 409 for drinking water (INN 2005) to evaluate the level of impact on water quality. Both NCh 1333 and NCh 409 consider total concentrations, which is consistent with the DGA data used in this work.

General descriptive analysis
Discharges in the Choapa River at Cuncumén (E4, upper zone) were the highest in the historical record (around 4 m 3 /s), which positions this river as the main source of surface water in the basin. The lowest values were observed in the Cuncumén River (E3) and in the upper part of the Chalinga River (E2), on the order of 0.2 and 0.5 m 3 /s, respectively, which also show the least variability. The most variable flows occurred in the middle and lower part of the Choapa River (E6, E8, and E9) ( Fig. 3(a)).
In general, the basin's rivers have maintained rather stable pH levels over time (i.e. ranging between 7 and 8.5), with most stations showing a slight increase in pH during the periods considered (see Supplementary material, Fig. S1). Also, whereas this parameter generally meets the water quality standards, it is possible to note a slight increase in pH from east to west throughout the basin (i.e. from higher to lower areas) (Fig. 4).
EC and major ion concentrations at all stations (Figs 3(b), 4 and Supplementary material, Fig. S2) decreased slightly for the period 1960-1990, and increased slightly afterwards. In general, the values tend to increase as one goes downstream in the basin, with the Illapel River (E8) having the highest (median) values. A possible explanation refers to the contributions of the Aucó Creek, where there is an important concentration of small mining operations (Parra et al. 2011). With respect to regulations, only Cl − , and SO 4 2− are included in both NCh 1333 and NCh 409; EC is in NCh 1333, whereas Mg 2+ is in NCh 409. Of these parameters, only in the case of EC and SO 4 2− did a few specific samples exceed the thresholds of NCh 1333 and NCh 1333 and 409, respectively.
Considering pH, EC and major ions simultaneously highlights the particular situation of the Cuncumén River station, E3 (Fig. 4). Indeed, this station presents the highest pH values (medians) among all the higher elevation (headwater) stations (E1 to E4). Also, and although the extension of the data series for this station is rather limited, EC and major ions (e.g. for Ca 2+ and SO 4 2− ) both show very high values when compared to the other headwater monitoring stations, and a different pattern than most of the stations, initially increasing and reaching a peak around 2004, when the values started to decrease ( Fig. 3 and Supplementary material, Figs. S3 -S8). The reason that the behaviour observed at E3 is very different from the rest of the monitoring stations is likely due to particular lithological characteristics, i.e. the upstream existence of moderate (i.e. propylitic) hydrothermal alteration zones (Parra et al. 2011), as well as mining activities associated to the MLP copper porphyry deposit, initiated in 1999. Parra et al. (2011) observed similar behaviour in SO 4 2− concentrations peaking in 2003-2004, which they explained by the combined effects of increased mining activity (MLP) and low precipitation in previous years. In addition, during this time, there was damage to rock glaciers Brenning 2008, Brenning and which are known to locally impact stream water quality (Rodriguez et al. 2016). The subsequent decrease in EC and major ion concentrations at E3 after 2004 could be the result of improvements in the company's environmental management, which was initiated in that year to gain approval to expand production capacity to 175 kt/d. These measures considered practices to avoid water contamination from tailings deposits (rock drains, settling ponds, etc.) and a commitment to constant monitoring of surface water quality in the watershed (CONAMA 2004). In addition, MLP has achieved a high degree of water reuse and recycling in response to its need to reconcile the scarcity of available water with its increased mining production.
Nutrients (N-NO 3 − in Fig. 3(c); P-PO 4 3− in Supplementary material, Fig. S9) generally showed greater variability than that observed for major ions. Some N-NO 3 − peaks occurred in 1992-1993, 1997-1998 and 2006-2007, which corresponded to wet years following periods of low rainfall. Also, the concentrations of N-NO 3 − and P-PO 4 3− are slightly higher in the middle and lower parts of the basin (i.e. E6 and E8), and in general there is a relatively similar variability throughout the basin. Of both constituents, only N-NO 3 − is included in the NCh 409, with a limit of 10 mg/L that was not exceeded during the period considered (despite the agricultural activity in the basin).
Regarding trace elements, the small number of stations with sufficient data (Table 2) presented some difficulty for quantitative comparisons. However, for the Choapa River at Cuncumén (E4), the only station with data for Zn, Al, Mn, and Fe, a general pattern is observed (Fig. 3(d), Supplementary material Figs. S10 -S15), with the highest values occurring by 1997 and 1999, after which the concentrations decrease only to increase and peak again in 2015 and 2017. The Fe concentrations exceeded the NCh 1333  and especially the NCh 409 on several occasions throughout the basin. Also, the upper part of the Choapa basin (E4 again) had the highest concentrations for the selected elements. In the case of Cu, higher values were observed at the Cuncumén River (E3), which is reasonable given the aforementioned existence of MLP and hydrothermal alteration zones in the area (Parra et al. 2011). Also, specifically regarding Fe, and in addition to E4, the Choapa River at Salamanca station (E6) presents high values as well (very similar to E4 in terms of both maximum recorded and variability). Finally, As concentrations tend to be high, especially in E4, frequently exceeding the threshold of the NCh 409.

C-Q relationships
As in the case of the time series data, scatter plots and regressions for selected parameters are presented ( Fig. 5; the rest of the cases are included in Supplementary material, * Slope significantly different from 0 (p value < .05); ** slope significantly different from 0 (p value < .01). Figs. S16 -S31). In addition, Table 4 presents the values of the regression slopes for each parameter at each station and their statistical significance, and Fig. 6 presents the "b" vs. CV C /CV Q plot.
On the basis of the regression slope ("b") values (Table 4) it is seen that overall, with the exception of E3, pH shows a constancy behaviour (i.e. with slopes that are not statistically different from 0), whereas EC and major ions show mostly a dilution behaviour. The possible reasons for the different behaviour of E3 (as discussed above, i.e. the effect of MLP), is consistent with Bonta (2005) who, in three small experimental basins in Ohio, USA, showed changes in C-Q relationships as a consequence of disturbance related to coal mining and reclamation practices. Rose et al. (2018) also indicated that the presence of exogenous sources or factors can locally alter C-Q relationships. In the case of nutrients (P-PO 4 3− , N-NO 3 − ), constancy behaviours are generally observed, which vary from chemodynamic to chemostatic, respectively. This is similar to results described by Musolff et al. (2015), Hunsaker and Johnson (2017) and Peña-Guerrero et al. (2020). Finally, in the case of trace elements, and despite limited data availability, enrichment and chemodynamic behaviours are predominantly identified for Fe, Mn, and Al, elements whose ionic potential makes them more susceptible to hydrolysis and which, therefore, tend to be incorporated in molecular form into the fluvial current. Indeed, for the pH range of the study area, it is normal that these constituents are not dissolved but are transported adsorbed to particulate material present in the water column. It is, then, expected to find higher total concentrations (consistent with the DGA data) associated with higher flows, which in turn are related to increased bottomsediment removal and suspended sediment loads being transported by the rivers (Cidu and Frau 2009). In contrast, Cu and Zn showed mixed enrichment and constancy behaviours (E3 and E4, respectively), with greater variability (i.e. chemodynamic pattern) for the former. Finally, As showed heterogeneous behaviour (i.e. dilution, constancy, enrichment).

Climate change effect
For the period 2010-2040, and based on the estimation of a decrease in Q of about 14%, the greatest variation in concentration (with respect to the average of the historical dataset) was found for trace elements, particularly Fe, Mn and Cu ( Fig. 7(a)). Based on the enrichment-type relationship identified above, we estimate a resulting reduction in C (with respect to the long-term average) on the order of 10% or lower. For parameters with a dilution-type behaviour, such as SO 4 2− , Cl − , Na + , Mg 2+ and B, their concentrations will likely increase, but again this change should not exceed 10%. For the period 2040-2070, considering a Q decrease by 20%, we estimated that concentrations of SO 4 2− and Cl − could increase to more than 10%, in particular in E4 ( Fig. 7(b)). Likewise, during this period, more stations are projected to show increases in major ions or EC by more than 5%, when compared to the 2010-2040 period. Of course, these changes ccould be greater if flow reductions are greater than forecasted values. Finally, it is interesting to note that E4 and E1, and to some extent E2 (i.e. headwater stations), will experience the greatest increases in major ions and salinity, since this Figure 6. CV C /CV Q vs. b relationships (average ± one standard deviation for all stations) for the different parameters considered (brown: field parameter; blue: major ion; green: nutrient; red: trace element).
could have a detrimental effect, especially on the agricultural activity carried out downstream.
When compared to the water quality standards, it is likely that average concentrations of Fe, As and Mn would still not comply with NCh 409, and eventually with NCh 1333 (e.g. in E4). Indeed, on one side, the enrichment behaviours identified for these parameters, especially Fe and Mn, suggest that lower Q (as expected from climate change) should translate into lower C. However, these lower concentrations will still be exceeding, on average, regulation thresholds (as inferred from the C-Q relationships based on historical data). Thus, special attention should still be given to the evolution of these constituents (and therefore related water quality monitoring) in the coming years, as well as to any given human action that may increase their presence in water (e.g. related to mining or other activities).

Addressing sources of uncertainties and their implications
As expected in hydrology-related studies, especially when considering future conditions (i.e. related to climate change predictions), there are sources of uncertainty worth considering for a robust assessment of the results. Thus, four sources of uncertainty are addressed, as follows: (a) For some parameters and stations, data availability was relatively low due to the quarterly/bi-monthy monitoring frequency by the DGA. Thus, it cannot be ruled out that intra-annual changes in flows may have translated into changes in the concentrations of certain constituents (Zhu et al. 2018, Montserrat et al. 2022. Indeed, given the scope of this study, the analysis presented corresponds to an expected average ("long-term") concentration of certain constituents in the rivers under future conditions. Certainly, in specific events or periods (e.g. years or months within a year), particular changes in the expected concentrations may differ from the levels derived from the C-Q relationships. (b) The time series for different stations differed across the basin. For example, stations E5 and E7 only had data to generate the C-Q regressions up to 1980 and 1977, respectively, while the rest of the stations had data up to 2019. This likely introduced a factor of uncertainty when projecting concentration variations for future flows, because one of the underlying assumptions of the analysis is that the C-Q relationship (type of behaviour) would remain relatively stable over time. (c) Although the slope "b" of the C-Q relationship (whether or not it was statistically different from 0), and the CV C /CV Q index are considered together as suitable behavioural classification criteria, we noted a high level of dispersion in the scatter plot for several parameters and stations, causing relatively low R 2 values for the C-Q relationships. Hunsaker and Johnson (2017) and Zhu et al. (2018) experienced similar issues with their datasets, which they related to other factors, besides Q, that may have influenced the constituent concentrations. In particular for trace elements, pH would be an important one. In fact, effects of climate change on water quality, especially in terms of trace elements, could be related to other factors beyond hydrological changes (streamflow). Although the climate projections considered in this work took into account the combined effect of changes in precipitation and temperature on Q, changes in temperature patterns may also affect various geochemical processes such as sulfide oxidation reactions in mineralized zones, the generation of acid drainage, and thus the solubility and transport of trace elements (Zarroca et al. 2021, Montserrat et al. 2022). (d) The existence of two streamflow regulation works (El Bato and Corrales reservoirs) may seem, at first glance, to represent an additional element of uncertainty in some of the C-Q relationships. El Bato's effect likely was minimal because it became operational only recently. Corrales' effect on Q likely was minor as well, because it is located in a tributary of the Choapa River (Camisas creek), and not along the main river course. Finally, even if there is some effect of the presence of these reservoirs, more than half of the monitoring stations in the headwaters (E1, E2, E3, E4, and E6) are clearly not affected by them.
Despite the issues described, the methodological approach yielded confident and consistent results about the general trends for the different constituent concentrations as a consequence of discharge reductions. Therefore, the approach in this study could be replicated in other watersheds, both in Chile and elsewhere, with similar challenges in determining likely water quality modifications due to climate change, in particular streamflow, even in similar conditions of data constraints.

Conclusions
We developed C-Q relationships of several parameters and constituents of rivers of the Choapa basin (Illapel, Chalinga and Cuncumén rivers), north-central Chile. Historical streamflow and constituent concentration data from nine monitoring stations operated by the DGA were utilized. Simple regressions were generated on the values of the C and Q data. Unlike most related work that is based exclusively on the characterization of C-Q relationships in either weathering products (or geogenic solutes) or trace elements, the present work considered a wide range of constituents. Also, and as a novelty of this work, we extended the use of the C-Q relationship to the assessment of possible effects of climate change on water quality. For most of the selected stations, we determined a constancy behaviour for parameters such as pH, N-NO 3 − , P-PO 4 3− , B, As, and Cu, although these last three presented dilution and enrichment patterns in some stations as well. The EC and major cations and anions showed dilution behaviour, with Na + , Mg 2+ , Cl − , SO 4 2− and K + showing the greatest variability with respect to streamflow variations. On the other hand, Fe, Al, Mn and Zn showed enrichment behaviours, with Al and Fe showing the greatest variability with respect to streamflow.
The behaviour of the parameters also varied with respect to the location of the stations. The data for the Cuncumén River (E3) showed a very different pattern from the other stations in most of the parameters. This is most likely due to the influence of local lithologies and particularly the presence of the MLP mining operation, which affected the behaviour of the different parameters from 2001 to 2004, as identified in the time series.
Regarding the projections of the effect on changes in streamflow, and subsequently, on the concentrations of the different constituents, the estimated percentages of variation of the latter did not exceed ±10% for the period 2010-2040, and ±15% for the period 2040-2070. Of the parameters that showed a dilution behaviour, SO 4 2− and Cl − were projected to have the greatest percentage increase. However, when applying these projections to historical averages, the expected concentrations were far from exceeding the limits of the water quality standards considered. In contrast, for parameters that showed enrichment behaviours, Fe, Mn, and Cu were estimated to have the greatest decrease in concentrations. However, given the historical concentrations, Fe and Mn remain water constituents that require special attention in monitoring efforts in the coming years.
Despite some sources of uncertainty because of data constraints, the proposed methodology is novel and practical for estimating the possible effect of climate change on the water composition in rivers such as those in the Choapa River basin, given forecasted changes in streamflow, provided it is possible to develop suitable C-Q relationships. In addition, the use of C-Q relationships and associated spatial and temporal variability (e.g. the same parameter between different monitoring stations or the same station between different periods) highlights it as a useful tool for the identification of particular effects of human activities or special environmental factors.