Reynald Eugénie

Averil Prost

Seth Taylor

Rita Zantout

On a workaround for an overflow in streaming process mining

This document serves as a report of a SEME project held in Pointe-à-Pitre, Guadeloupe, during June 2023.

Introduction

Context

It is a sacred rule in programming that each action must leave a trace. Computers and communication devices are constantly recording their states and various failures in log files, with a high level of detail. These flows of information stay widely unexploited, due to their technical complexity and overwhelming volume. However, they contain invaluable a posteriori observations of the interactions of applications, in particular regarding the eventual flaws or malfunctions.

The intelligent treatment of log files is the starting point of process mining. This recent field is developed both in the academic community and in the industry. In the latter case, the target technique is the streaming process mining, in which the generated data is immediately handled. This is the main service offered by Softbridge Technology, through a platform implemented in Java.

One of the aims of Softbridge is to detect latency in a process. The proposed method implies to approximate the distribution of durations of the events contributing to this process, such that outliers may be identified. However, the choice of an on-line treatment of incoming data puts a constraint on the representation of these durations, since each event must be processed only once. Consequently, some basic-but-efficient methods such as histogram representations are illadapted, since the range of datas is not known a priori. Instead, one can use a specific class of parametric representations, under the assumption that the parameters may be updated on-line. It is the case, for instance, of the generic estimators of the mean and standard deviation of a given probability law (see Section 1.2).

This project arises from a technical difficulty encountered in the implementation of the parametric approach. The duration of a process, denoted x i , is expressed in milliseconds (following the classical POSIX representation). Some processes may run up to 10 months, and the associated x i becomes so large that x 2 i overflows the memory allocated by the primitive type long of Java. The language implements a class BigInteger that virtually achieves unbounded precision, but at the expense of a greater computational cost and possible loss of portability. The purpose of this document is to find the most efficient workaround, respecting the constraints of streaming estimation, using only primitive types, and maximizing the precision.

Formulation of the problem

In the sequel, the processes will be observed only through their duration in milliseconds, generically denoted x i . Suppose that N observations are available, and denote

µ N = 1 N N i=1 x i , S N = 1 N N i=1 x 2 i ,
as the mean and mean of the squares respectively. Our aim is to parametrize the distribution of processes by either a normal or a log-normal law. Recall that the normal and log-normal distributions are given by

f [µ, σ](x) = 1 √ 2πσ exp - (x -µ) 2 2σ 2 , and g[µ, σ](x) = 1 x √ 2πσ exp - (ln(x) -µ) 2 2σ 2 ,
An on-line estimation of the relevant moments to compute and plot the distribution is given by the following algorithm:

Algorithm 1: Naive algorithm

1 Set N = 0, µ 0 = 0 and S 0 = 0. 2 while the program is running do 3 Get a new duration x N +1 .
4

Update the mean by µ N +1 := (N × µ N + x N +1 ) /(N + 1).

5

Update the sum of squares by

S N +1 := N × S N + x 2 N +1 /(N + 1). 6 if the law is evaluated then 7 Compute σ N +1 = S N +1 -µ 2 N +1 .
8

Compute and plot the chosen distribution. 9 end 10 end

In practice, N • S N is the value stored and input into the formula for the standard deviation σ. This poses two issues associated with overflow in this algorithm, both occurring at line 5.

Problem 1. If the process has duration x i > 2 32 -1 then overflow will occur when squaring.

Problem 2. The sum of squares N • S N can become unbounded and will cause overflow if stored directly.

Proposed solutions

In theory, Java has a non-primitive type called BigInteger which allows for an adaptivity from the int type into a long type when an overflow is detected. In practice however, this solution is not satisfactory since this type is not generic and its storage in memory is not necessarily controlled. Furthermore, Softbridge seeks a solution methodology which is independent of the language and equipment. We thus propose two solutions which address the problems posed by Problems 1 and 2 which use only integer or floating-point encodings. Each method is explained and detailed in the following sections. We then characterize the possible advantages and drawbacks realized using these approaches.

Emulation of the behavior of BigInteger

The first proposed approach consists in emulating the operations of the BigInteger type using only operations on the int type until the variance is accessed. The central idea behind this approach is to use modular arithmetic to have a larger representation of an integer greater than 2 63 -1. We consider an emulation of the behavior of the BigInteger type of Java, by decomposing the duration using two values A, B, each stored as primitive types.

Following the Java convention of signed natural numbers, we define the following sets.

Definition 1 (Int64). A natural number x admits a representation as an Int64 if it belongs to -2 63 , 2 63 -1 .

Definition 2 (MyBigInt). A natural number x admits a representation as a MyBigInt if it is written as x = A × 2 e + B, where A, B are Int64, and e ∈ 6, 24 .

The range defining the exponent will be clarified based on the operation carried out using this solution. Note that, the set of MyBigInt does not coincide with the set of numbers represented in Java as BigInteger.

Assumption 1 (Range of incoming data). In the following, we will assume that the incoming values (x i ) i are all exactly represented as Int64.

As a consequence of Assumption 1, the mean values (µ N ) N are all exactly represented in 64-bits integers.

Assumption 2 (Range of the squares). We will also assume that at any given N , the sum N i=1 x 2 i has a representation as a MyBigInt.

Let // denote the integer division, i.e. a//b = ⌊a/b⌋, and mod be the modulo operation (for completeness, let us point that we consider the multiplication to have priority upon the modulo, and denote ab mod c = (ab) mod c).

Remark 1 (Range of validity of the assumptions). Assumptions 1 and 2 are based on the following computations: a process lasting 4 years will take α = 4 × 365 × 24 × 3600 × 1000 = 126144000000 ≃ 2 36 ms to complete. This number still admits a representation as an 64-bits integer. On the other hand, the achievable precision of a MyBigInt is of the order of magnitude of (2 2e -1)2 2e ≃ 2 4e . Hence, in the worst case where all the processes last 4 years and taking e = 24, the number of processes N such that N × α overflows the MyBigInt representation is at least

2 4×e //α 2 + 1 = 4 979 048.
In comparison, the number of 4-year processes needed for the sum of squares to overflow the double representation is around 10 285 . We consider this to be a reasonable amount.

The manipulation of MyBigInt relies on the following lemma.

Lemma 1 (Operations). Let x ∈ R + be given by x = a2 e + b, where a and b belong to 0, 2 e -1 . Then its square x 2 may be written as A × 2 2e + B, where

A = a 2 + 2ab2 e + b 2 //2 2e , and B = (2ab2 e mod 2 2e + b 2 ) mod 2 2e .
(1)

Let now x i = A i 2 2e + B i , with A i , B i ∈ 0, 2 2e -1 and i ∈ {0, 1} .
1) The sum x 0 + x 1 writes A2 2e + B, where

A = A 0 + A 1 + (B 0 + B 1 )//2 2e and B = (B 0 + B 1 ) mod 2 2e .
2) Assume that x 0 ⩾ x 1 , and denote s = 0 if B 0 ⩾ B 1 , and s = 1 if B 0 < B 1 . The difference x 0 -x 1 writes as A2 2e + B, where A = A 0 -A 1 -s, and B = B 0 -B 1 + s 2 2e .

3) Let n ∈ N such that nB 0 can be represented as an Int64. Then y = nx 0 may be represented as C2 2e + D, where C = nA 0 + (nB 0 )//2 2e , and D = (nB 0 ) mod 2 2e .

Remark 2 (Choice of the exponent). When calculating 2a i b i 2 e , there is a risk of overflowing. We can circumvent this issue by considering an adapted value for the exponent e. If we assume that any given duration is less than four year (∼ 2 37 ) then we have that a ∈ [0, 2 37-e ] =⇒ a 2 ∈ [0, 2 74-2e ]. Since we require 74 -2e < 63, as a lower bound we have e > 5. As an upper bound, we note that since a ∈ [0, 2 37-e ] and b

∈ [0, 2 e -1] then 2 • a • b • 2 e ∈ [0, 2 38+e
] hence we require e < 25. In practice we take e ∈ [6, 24].

Proof. Notice that

x 2 = (a2 e + b) 2 = a 2 2 2e + 2ab2 e + b 2 .
Recall that ab mod c = ((a mod c)(b mod c)) mod c. The expression (1) follows from

B = x 2 mod 2 2e = 2ab2 e + b 2 mod 2 2e = (2ab2 e ) mod 2 2e + b 2 mod 2 2e ,
and, since all the terms are nonnegative,

A = x 2 //2 2e = a 2 2 2e + 2ab2 e + b 2 2 2e = a 2 + 2ab2 e + b 2 //2 2e .
Item 1) is straightforward. Item 2) may be treated by disjunction, and Item 3) follows the same reasoning as Item 1).

To completely re-implement the naive algorithm Algorithm 1, one needs to have an approximation of the square root of a natural number x = A2 2e + B. Our implementation amounts to compute

√ x ≃ A + B/2 2e × 2 e , (2)
where the first part is encoded as a double. A lot of submethods may arise from this point, by considering the order of operations and the point to which the encoding falls back into double. The current implementation turns B into a double, then performs the division B/2 2e in floating-point, and converts A to a double on the way.

Remark 3 (Alternative). Historically, we considered an approximation as A 0 * 2 2e /(A 0 + 1) + B 0 /(A 0 + 1) × √ A 0 + 1 in order to reduce the magnitude of the stored numbers, but the results turn out to be unsatisfactory (see Section 3.2).

Altogether, our first solution strategy can be implemented as the following pseudocode algorithm:

Re-scaling using floating-point encoding

As a second proposed solution, as initially hinted by Softbridge, we consider a dynamic re-scaling of the sum of squares.

Let λ > 0 denote our re-scaling parameter and let σ 2 λ := λσ 2 . A re-scaled distribution takes the equivalent form

f [µ, σ](x) = 1 √ 2πσ exp - (x -µ) 2 2σ 2 = √ λ √ 2πσ λ exp -λ (x -µ) 2 2σ 2 λ =: f λ [µ, σ λ ](x).
In the above expression, λ may be dynamically chosen small in order to keep the value σ λ in ranges of values that may be represented in floating-point encoding (see 4) with sufficiently high precision. Note here that the terms λx 2 i have no reason to stay natural numbers and that by taking λ = 1 we recover the case where each number is simply turned into a float.

Remark 4 (Double-Precision). The IEEE standard floating-point encoding of double precision uses 64 bits, broken down into the following sections:

• a first sign bit s.

• 11 bits of exponent, encoding a number p ∈ 0, 2 11 -1 .

• 52 bits of mantissa, that we will label in decreasing order as b

51 , b 50 , b 49 , • • • , b 0 , encoding a multiplier m ∈ [1, 2[ by m = 1 + 52 i=1 b 52-i 2 -i .
The floating-point double will then be represented as

x = (-1) s × m × 2 p-1023 = (-1) s × 1 + 52 i=1 b 52-i 2 -i × 2 p-1023 .
We will say that a real number x is encodable if it may exactly be represented in double floating-point encoding. Given a fixed exponent p, the encodable values are equally distributed in the interval [-(2 -2 -52 )2 p-1023 , (2 -2 -52 )2 p-1023 ], with a step of 2 -52 × 2 p-1023 .

Re-scaled operations

The computation of the mean remains unchanged, however the computation of the sum of squares becomes

S λ N +1 = N × S λ N + λ x 2 N +1 N + 1 = N N + 1 × S λ N + 1 N + 1 √ λ x N +1 2 .
The computation of the re-scaled deviation at step N will be

σ 2 λ = S λ N -( √ λ × µ N ) 2 .
The last modification of the naive algorithm will intervene in the evaluation of the law, that has to take into account the parameter λ. The modification will depend on the particular expression of the distribution. The correct expression of the re-scaled normal and log-normal laws are the following:

f λ [µ, σ λ ](x) := √ λ σ λ √ 2π exp -λ (x -µ) 2 2σ 2 λ and g λ [µ, σ λ ](x) := √ λ xσ λ √ 2π exp -λ (ln(x) -µ) 2 2σ 2 λ . (3)
Implementation as an on-line algorithm Given the wide range of possible durations, the choice of λ is somewhat arbitrary. We consider two possible implementations:

• Fixing λ a priori, allowing us to detect whether a small or a large value of λ is well-suited.

• Update λ dynamically as the inverse of the current maximum of x i , or the current maximum of x 2 i , reducing the order of magnitude of the data.

If we denote λ N as the choice of λ at step N , an algorithm for the proposed second solution strategy is as follows:

Algorithm 3: Re-scaling algorithm

1 Set N = 0, µ 0 = 0, S 0 = 0 and λ 0 = 1. 2 while the program is running do 3 Get a new duration x N +1 . 4 Update the mean by µ N +1 := N N +1 × µ N + x N +1 /(N + 1).
5

Update the scaling parameter λ N +1 (while still storing λ N ).

6

Update the sum of squares by S

λ N +1 N +1 := S λ N N × λ N +1 λ N + λ N +1 x 2 N +1 . 7 if the law is evaluated then 8 Compute σ N +1 = S λ N +1 N +1 /(N + 1) -λ N +1 µ 2 N +1 .
9

Compute and plot the re-scaled chosen distribution as in (3). 3 Numerical results

Testing environment

In order to assess the different proposed solutions, we designed a range of test cases. We divide them into deterministic test cases, in which every parameter is analytically known at every time, at the expense of being non-realistic, and random test cases that mimic the industrial context. A test case is determined by 1. a method to provide the sequence (x i ) i∈N , 2. a set of expected parameters µ N and σ N in function of the number N of available processes.

Let μN (resp. σN ) be the numerical approximation of µ N (resp. σ N ). As a measure of the error for the proposed algorithms, we consider the following indicators:

• a measure of the overflow, equal to 0 if the program is completed without overflowing, and +∞ otherwise,

• the relative error in mean e µ := |μ N -µ N | /µ N ,

• the relative error in standard deviation e

σ := |σ N -σ N | /σ N ,
The reported results do not show problems of overflow, and we do not mention it.

Deterministic test cases

As deterministic test cases, we consider the following two trials 1. the constant test case: x i = x a fixed value. The expected parameters are µ N = x and σ N = 0 at any step N .

2. the cyclical test case: the values of x i will cycle through a given finite set X = {x 1 , x 2 , • • • , x m }. Let µ be the mean of X, and S = 1 m m i=1 x 2 i . The expected parameters are given by

µ N = 1 N N i=1 x i = N//m N mµ + 1 N N %m i=1 x i , and σ N = 1 N N i=1 x 2 i -µ 2 N = N//m N mS + 1 N N %m i=1 x 2 i -µ 2 N .

Random test cases

In the random case, the samples x i are drawn from a given law. Consequently, the values of µ N and σ N will only be approximated by their numerical for a large number of samples. We keep the law fixed with respect to N , so that the expected parameters µ = µ N and σ = σ N are independent of N . We consider the following two trials.

1. The uniform discrete distribution U(1, K) between 1 ms and K ms, where K is even and is taken as a large duration of ≈ 10 months. The expected parameters are given by µ = (1 + K)/2, and the standard deviation is σ = k 2 -1 12 .

2. The normal distribution N (µ, σ) for known mean and standard deviation.

Comparison of methods

In all the following, the number of simulated processes is taken as N = 100 for deterministic test cases, and N = 10 000 for random ones. We compare the different methods against an implementation using BigInteger. The plotted errors are computed as the mean of relative errors over the last 10 iterations. We begin by some very simple deterministic cases to assess the method. The first test case (represented in Figure 1) considers a constant sequence of durations x i ≡ x. The value of x is set to 2 i , with i varying in 8, 38 in order to make the sum of squares S overflow the long representation. The x axis represents the natural number i ∈ 8, 28 such that x = 2 i . The y axis shows the relative error on, respectively, the sum of the squares S = N i=1 x 2 i , and the standard deviation σ = S/N -µ 2 . The method scaling max corresponds to the choice of λ N = max(λ N -1 , x N ), with λ 0 = 0. The method float corresponds to the constant choice of λ = 1. The method emul 1 corresponds to the Algorithm 2, with the approximation of the square root performed as in (2). Unrepresented methods are giving far larger errors (see Tables 1 and2).

The second test case (see Figure 2) considers a sequence x i alternating between 1 and the value 2 i . The i is chosen increasingly as before.

To conclude with analytically known test cases, the Figure 3 considers a periodic sequence over a fixed set of values, whose maximum is again set to 2 i , with i varying as before. The series of data labeled scaling min is an instance of the scaling algorithm for a fixed value λ = 1e -12. The series emul 2 refers to the emulation of BigInteger behavior, with an implementation of the square root with a division by the largest factor (see Remark 3). The remaining labels are identical to Figure 1.

To get closer to reality, we also considered random test cases. The error represented in the following figures is computed relatively to the mean and standard deviation of the laws that we sample, and the reader may keep in mind that part of the error may be due to the random sampling. Figure 4 presents the of the uniform distribution in [1, 2 i ], with i increasing as previously.

Figure 5 shows the sampling of a normal law of increasing mean and fixed standard deviation at 200. This case is the most relevant in the given context.

At first sight, it seems natural that the scaling method should give better results for a suitably adapted parameter λ, that would keep the stored quantities into reasonable intervals. To assess this, Figure 6 represents the variation of the relative error of S and σ for a given λ, ranging from 1 to 10 -38 . The test case is chosen to be the alternation between 1 and 2 8 , so that the error is not influenced by a random sampling.

The comparison is summarized with better precision in the following tables, in which the numbering of simulation coincide with that of the corresponding figures. For each x ∈ 8, 37 , the values are taken periodically in the set {1000, 2 x , 6000, 3000, 5 × 2 x//2 }. The scaling parameter λ is chosen as 10 -x , where x ∈ 0, -38 lies on the first axis.

Conclusion

Our project explored two methods to avoid overflowing in a stream process mining used by Softbridge technologies, namely 1. emulating BigInteger with long types: decompose a large number in order to store it as A × 2 e + B, where A and B may be represented as long values, 2. re-scaling the parameters: store a scaled version of the parameters, thus converting any long to double at the beginning of the process.

In regard of the first method, we observed that we could represent the sum of the square without any loss of information. However, for the computation of σ, there is a non-negligible loss of precision associated with the square root operation. We have tested many options to calculate the square root, and presented the results for some factorisations. Improving the computation of the square root may allow more precise results of the σ.

On the second method, the loss of precision occurs in the computation of the sum of the squares. Nevertheless, this loss is mitigated when calculating the standard deviation, showing better results than method 1. We also observed that the choice of the scaling parameter λ does not have a significant influence on the results.

In conclusion, our advice to the Softbridge company would be to turn to a double encoding. This choice is by far the most easily implemented, and the above study shows that it may outperform more refined methods.

Algorithm 2 :69

 2 Emulation of BigInteger algorithm 1 Set N = 0, µ 0 = 0 and S 0 = 0. 2 while the program is running do 3 Get a new duration x N +1 . 4 Update the mean by µ N +1 := (N × µ N + x N +1 ) /(N + 1). 5 Decompose x 2 N +1 in MyBigInt representation. Perform the sum S N +1 := S N + x 2 N +1 in MyBigInt representation using the formula (1). 7 if the law is evaluated then 8 Decompose (N + 1) × µ 2 N +1 in the MyBigInt representation (using a truncation). Perform the difference A2 2e + B := S N +1 -(N + 1)µ 2 N +1 in MyBigInt.

  to obtain σ.16Compute and plot the chosen distribution.

  17

Figure 1 :

 1 Figure 1: Relative error in S (a) and σ (b), constant test case.

Figure 2 :

 2 Figure 2: Relative error in S (a) and σ (b), alternating in {1, 2 i }.

Figure 3 :

 3 Figure 3: Relative error in S (a) and σ (b) alternating in a broader range of values.

Figure 4 :

 4 Figure 4: Relative error in S (a) and σ (b), uniform distribution in [1, 2 i ].

Figure 5 :

 5 Figure 5: Relative error in S (a) and σ (b), normal distribution with mean 2 i and standard deviation 200.

Figure 6 :

 6 Figure 6: Relative error in S (a) and σ (b), alternating test case and varying λ.

Table 1 :

 1 Errors for the sum of squares S

	Simulation n°error scaling max	float	scaling min	emul 1	emul 2
	1	max	3.88e-17	3.88e-17	4.61e-16	0.00e+00 0.00e+00
	1	mean	1.12e-17	1.12e-17	4.30e-16	0.00e+00 0.00e+00
	2	max	1.19e-15	1.19e-15	2.20e-15	0.00e+00 0.00e+00
	2	mean	8.33e-17	8.33e-17	1.05e-15	0.00e+00 0.00e+00
	3	max	1.03e-15	1.03e-15	9.70e-16	0.00e+00 0.00e+00
	3	mean	1.06e-16	1.06e-16	3.87e-16	0.00e+00 0.00e+00
	4	max	2.24e-02	2.24e-02	2.24e-02	0.00e+00 0.00e+00
	4	mean	6.98e-03	6.98e-03	6.98e-03	0.00e+00 0.00e+00
	5	max	7.64e-03	7.64e-03	7.64e-03	0.00e+00 0.00e+00
	5	mean	6.20e-04	6.20e-04	6.20e-04	0.00e+00 0.00e+00

Table 2 :

 2 Errors for the standard deviation σ

	Simulation n°error scaling max	float	scaling min	emul 1	emul 2
	1	max	0.00e+00	0.00e+00	2.90e+03	0.00e+00 0.00e+00
	1	mean	0.00e+00	0.00e+00	1.93e+02	0.00e+00 0.00e+00
	2	max	2.17e-07	2.17e-07	2.17e-07	2.61e-03	nan
	2	mean	2.17e-08	2.17e-08	2.17e-08	2.25e-04	nan
	3	max	9.21e-08	9.21e-08	9.21e-08	3.41e-04	nan
	3	mean	5.26e-09	5.26e-09	5.26e-09	5.70e-05	nan
	4	max	9.00e-03	9.00e-03	9.00e-03	1.35e-02	nan
	4	mean	2.85e-03	2.85e-03	2.85e-03	3.56e-03	nan
	5	max	4.75e+01	4.73e+01	5.49e+01	1.16e+03 1.15e+03
	5	mean	2.82e+00	2.82e+00	3.14e+00	1.82e+02 1.82e+02


Pauline Louis

Davide A L Vignati

Steve Pontvianne

Marie-Noëlle Pons
email: marie-noelle.pons@univ-lorraine.fr

Spatial distribution of Rare Earth Elements in a transnational watershed: the case of the Danube River

Keywords: Anomaly, Gadolinium, Magnetic resonance imaging, Ytterbium

The distribution of Rare Earth Elements was investigated along the Danube River (from its sources in Germany to its delta in Romania/Ukraine) and in some of its tributaries for the first time during the fourth Joint Danube Survey (JDS4) in 2019. The negative Ce anomaly and the La/Lu ratio tend to decrease between Ulm (Germany) and the Black Sea, indicating an increasing trend of the anthropogenic influence as the Danube River runs toward the sea. This diffuse influence is also attested by the slow spatial variation of classical pollution parameters such as dissolved organic carbon (amount and optical characteristics) and nitrates. A strong gadolinium (Gd) anomaly (between 6 and 10), due to the release of anthropogenic Gd, was detected starting 10 km downstream the sources of the Breg River and the Brigach River, whose junction forms the Danube River. The anomaly remains larger than 10 until the Danube River leaves Croatia at Ilok. The anomaly pattern is related to the difference in the use of magnetic resonance imaging (MRI) in the countries sharing the Danube River watershed. The estimated annual flux of anthropogenic gadolinium in the Black Sea is 0.9 ton/year. An Ytterbium anomaly was detected in stations close to the delta and needs further investigation.

Introduction

Rare Earth Elements (REEs), i.e. the fifteen lanthanides (from Lanthanum to Lutetium) to which Yttrium and Scandium are often associated, are increasingly present in our modern life: Y, La, Ce, Eu, Gd, and Tb in fluorescent lamps [START_REF] Hobohm | Optimized elemental analysis of fluorescence lamp shredder waste[END_REF], Ce in batteries [START_REF] Porvali | REE(III) recovery from spent NiMH batteries as REE double sulfates and their simultaneous hydrolysis and wet-oxidation[END_REF], Nd, Pr and Dy in permanent magnets, La in catalysts (Kulaksız and Bau, 2011a), Ce in diesel fuel additives [START_REF] Dale | Transformation of cerium oxide nanoparticles from a diesel fuel additive during combustion in a diesel engine[END_REF] are some examples. Anthropogenic uses of REE can both increase their concentrations above natural background levels and alter their natural distribution pattern by introducing anomalies in REE profiles. While the solid waste that may result from their use can be collected and potentially reprocessed [START_REF] Hobohm | Optimized elemental analysis of fluorescence lamp shredder waste[END_REF][START_REF] Yang | REE recovery from end-of-life NdFeB permanent mgnet scrap: a critical review[END_REF][START_REF] Porvali | REE(III) recovery from spent NiMH batteries as REE double sulfates and their simultaneous hydrolysis and wet-oxidation[END_REF][START_REF] Becci | Sustainable strategies for the exploitation of end-oflife permanent magnets[END_REF], the same is not true for REEs used in the medical field, in particular in contrast agents. These include Gd-based contrast agents (GBCAs), used in Magnetic Resonance Imaging (MRI) examinations. GBCAs are eliminated through the urine and are therefore found in urban wastewater. These complexes are nonbiodegradable and pass through wastewater treatment facilities to be released into aquatic environments. Since the first evidence (Bau and Dulski, 1996), the presence of Gd has been detected in aquatic ecosystems in many countries, even when they do not have MRI facilities.

Table S1 provides data extracted from literature about Gd concentrations in non-saline surface waters around the world. It is difficult to calculate an average value as the data were obtained over 25 years in different hydrological conditions, in pristine waters as well as in highly urbanized areas. The highest value was measured (492 ng/L) in Postdam (Berlin/Brandenburg Metropolitan Area) by Kulaksız and Bau (2011b).

Although Gd is probably the most documented REE as a micropollutant of increasing concern, anthropogenic La has been detected in the Rhine River in Germany (Kulaksız and Bau, 2011a).

In water samples (Kulaksız and Bau, 2013) as well as in shells [START_REF] Valdés-Vilchis | Rare Earth Elements and Yttrium (REE+Y) patterns in recent Anadara brasiliana shells from Playa Norte, Barra de Cazones (Veracruz, Mexico): Evidence of anthropogenic contamination linked to river output?[END_REF], samarium positive anomalies have been related to petrochemical industries, where Sm-based catalysts are used. Furthermore, the industrial activities related to the extraction and purification of REEs as well as REE-based products manufacturing can release REEs in aquatic ecosystems. According to [START_REF] Chen | The profiling of elements and pesticides in surface water in Nanjing, China with global comparisons[END_REF], contamination by REEs of the Qinhuai River in China might be related to electronic waste recycling and permanent magnet industries. In Dingnan County (China) heavy contamination of soil and water has been recorded close to a REEs mine [START_REF] Liu | Water, sediment and agricultural soil contamination from an ion adsorption rare earth mining area[END_REF]. Mining activities not focusing directly to REEs can also be the source of contamination by these elements, as shown by [START_REF] Khan | Assessing anthropogenic levels, speciation, and potential mobility of rare earth elements (REEs) in ex-tin mining area[END_REF] in a former tin mining area. Runoff from agricultural fields where municipal sludge [START_REF] Kaegi | Quantification of individual Rare Earth Elements from industrial sources in sewage sludge[END_REF] and/or chemical fertilizers [START_REF] Hu | Physiological and biochemical effects of rare earth elements on plants and their agricultural significance: a review[END_REF][START_REF] Turra | Sustainability of rare earth elements chain: from production to food -a review[END_REF] have been applied can participate also to surface waters pollution by REEs. Urban runoff, contaminated by particulate pollution related to coal burning and transport activities, is also a contributor of La and Ce discharges in surface water [START_REF] Shajib | Rare earth elements in surface specific urban runoff in Northern Beijing[END_REF].

Once anthropogenic REEs are detected in surface water, there is a risk of transfer to drinking water. Bank filtration is used in many countries for drinking water production and its efficiency as a barrier against transfer of pharmaceutical residues has been questioned [START_REF] Kondor | Occurrence of pharmaceuticals in the Danube and drinking water wells: Efficiency of riverbank filtration[END_REF]. As early as 2009, Gd has been detected as a micropollutant in tap water in Berlin (Kulaksız and Bau, 2011b), with a rapid increase over a few years [START_REF] Tepe | High-technology metals as emerging contaminants: Strong increase of anthropogenic gadolinium levels in tap water of Berlin, Germany, from 2009 to 2012[END_REF]. Bank filtration is not the only way to get anthropogenic Gd into tap water: karst aquifers can be easily contaminated as shown by Rabiet et al. (2009) in the Hérault watershed (France) and [START_REF] Boester | Utilize gadolinium as environmental tracer for surface watergroundwater interaction in Karst[END_REF] in Bavaria. Recently anthropogenic Gd has been reported in four tap waters in China, in the district of Guiyang (Han et al., 2021). If Gd-related health risks were thought to be minimal when GBCAs started to be used in the 1980s, more recent data indicate possible retention in different human organs such as kidney [START_REF] Fraum | Gadolinium-based contrast agents: a comprehensive risk assessment[END_REF], brain [START_REF] Chazot | Brain MRIs makeup the bulk of the gadolinium footprint in medical imaging[END_REF] or bones [START_REF] Veiga | Presence of other rare earth metals in gadolinium-based contrast agents[END_REF]. [START_REF] Thomsen | Are the increasing amounts of gadolinium in surface and tap water dangerous?[END_REF] has suggested monitoring more closely Gd as risks associated with tap water are not well known. Accumulation of REE has been detected in many aquatic organisms (such as crabs [START_REF] Lavezzo | Detection of rare-earth elements using fiddler crabs Leptuca leptodactyla (Crustacea: Ocypodidae) as bioindicators in mangroves on the coast of São Paulo, Brazil[END_REF] and bivalves shells [START_REF] Merschel | Rare earth elements in the aragonitic shell of freshwater mussel Corbicula fluminea and the bioavailability of anthropogenic lanthanum, samarium and gadolinium in river water[END_REF][START_REF] Goff | Compound-specific recording of gadolinium pollution in coastal waters by great scallops[END_REF][START_REF] Mouchi | Rare earth elements in oyster shells: provenance discrimination and potential vital effects[END_REF][START_REF] Valdés-Vilchis | Rare Earth Elements and Yttrium (REE+Y) patterns in recent Anadara brasiliana shells from Playa Norte, Barra de Cazones (Veracruz, Mexico): Evidence of anthropogenic contamination linked to river output?[END_REF] with the Gd anomaly linked to anthropogenic uses remaining detectable in bivalves [START_REF] Pereto | Corbicula fluminea: A sentinel species for urban Rare Earth Element origin[END_REF]) and fish [START_REF] Lortholarie | Rare earth element bioaccumulation in the yellow and silver European eel (Anguilla anguilla): A case study in the Loire estuary (France)[END_REF]. However, REEs seem to undergo biodilution along the food chain [START_REF] Amyot | Fate and trophic transfer of rare earth elements in temperate lake food webs[END_REF][START_REF] Santos | Fate and trophic transfer of rare earth elements in a tropical estuarine food web[END_REF] and the possible risks associated to their dietary uptake remain to be elucidated. Recently, [START_REF] Souza | Evidence of high bioaccessibility of gadolinium-contrast agents in natural waters after human oral uptake[END_REF] discussed the evidence of high bioaccessibility of gadolinium-contrast agents after ingestion of Gd-contaminated tap water.

If anthropogenic REEs have been detected in different aquatic ecosystems, in particular in rivers, the number of samples collected is often limited to a couple as it can be seen in Table S1 (Supplementary Information section). It has to be recognized that sampling in watersheds of large rivers requires sophisticated logistics: Zhang et al. (1998) could only investigate the Yangtze River lower section (eight stations on about 1,000 km of a 6,300 km-long river) downstream of the city of Wuhan. Ma and Wang (2023) managed to sample 15 stations on tributaries of the complex Pearl River delta. In such cases, it is difficult to relate the REE concentration to the watershed characteristics and disentangle the actual contributions due to anthropogenic pressures from geological factors (Louis et al., 2020).

The fourth Joint Danube Survey (JDS4), which took place in summer 2019, was an extraordinary opportunity to monitor for the first time the REEs along the Danube River and in some of its tributaries. The Danube River watershed is the second largest in Europe with a surface of about 800,000 km 2 and is shared by nineteen countries, making it the most international river basin in the world, with a total population of 81 million inhabitants.

Organized by the ICPDR (International Commission for the Protection of the Danube River)

Monitoring and Assessment Expert Group, JDS4 made possible the sampling in twelve countries over a period of about one month (JDS4, 2021).

The purpose of the present study is to report for the first time the spatial distribution of REE concentrations in the Danube watershed and put this distribution in relation with those of other elements (alkaline earth metals, major ions (i.e., sulfates, chlorides and nitrates)), dissolved organic carbon (DOC) concentration and typology and the population living in this transnational watershed.

Materials and methods

Sampling

The full list of all sample locations, with their codes, is given in the Supplementary Information section (Table S2 for the Danube River and Table S3 for its tributaries). Sampling upstream from Ulm took place July 3 rd (between Neu Ulm and Donaueschingen) and July 4 th 2019 (along the Breg River and the Brigach River between Donaueschingen and the river sources) (Figure 1). Samples (n = 35) were grabbed from bridges (except at the sources), 5 cm below the surface using a bailer made of high-density polyethylene (PEHD) (SDEC, Reignac-sur-Indre, France).

Two tributaries were also sampled (Eisenbach River and Iller River) and two drinking water samples were collected in Donaueschingen and Ulm. The samples were transferred into clean 250 mL HDPE (High-Density PolyEthylene) bottles and kept in a coolbox in the dark until the arrival in the lab. They were then stored at 4°C in the dark pending treatment within 24hrs. Sampling downstream from Neu Ulm took place in mid-July, 2019, within the framework of the JDS4 (Figure 2). The sampling team in each country had previously received clean 500 mL HDPE bottles. Fifty-one samples were grabbed along the Danube and some of its tributaries from the banks or from the middle of the river. Groundwater samples were also collected in wells recharged with Danube water through bank filtration in Austria, Slovakia, Hungary, Croatia, Serbia, Romania and Bulgaria (Figure 2).

All samples were stored in a coolbox in the dark and transferred to the Environmental Institute (Kos, Slovakia) prior to their overnight shipment to France where they underwent pretreatment within one day. and some of its tributaries (orange dots) and groundwater samples (yellow diamonds).

Analytical methods

Sample filtration and ancillary parameters

Samples were filtered through regenerated cellulose syringe filters (Phenomenex®) with 0.45 µm nominal pore size. Only the filterable fraction, which consists of dissolved ions and colloids < 0.45 μm, is considered in the present paper. Aliquots (50 mL for REEs, 10 mL for major elements) of the filterable fraction were acidified immediately after filtration with ultrapure HNO3 (at 1% v/v, acid class "Optima" for REEs and at 2% v/v, acid class "Trace Metal" for major cations). Other sample aliquots (80 mL) were preserved at 4 °C without acidification for measurements of dissolved organic and inorganic carbon (DOC and DIC), total nitrogen (TN), major anions and characterization of the dissolved organic matter (DOM) by optical methods (section 2.2.3). DOC and DIC were measured with a Shimadzu TOC VCHS combined with a TNM-1 module for TN determination. Major anions (Cl -, SO4 2-, NO3 -) were determined by ion chromatography (Thermo Scientific Dionex iCS 3000) and major cations (Mg, Na, K, and Ca) by inductively coupled plasma atomic emission (ICP AES) (Thermo Scientific iCAP 6300 Duo).

REEs analysis

REEs were determined using Inductively Coupled Plasma Mass Spectrometry (ICP-MS) (Thermo Scientific iCAPQ+prepFast) without preconcentration. Re and Rh at 50 ppb were used as online internal standards. The instrument was tuned for oxide and doubly charged ion formation, and interference corrections for oxide and hydroxide were applied when necessary. The analytical error was below 5%. The quantification limit, calculated as 10 times the standard deviation of replicated blank measurements, was 1 ng/L for all REEs. The filtration blank results (n = 15) were always below detection limits. SLRS-6 reference water was used to control the ICP-MS accuracy and reproducibility [START_REF] Yeghicheyan | A new interlaboratory characterisation of silicon, rare earth elements and twenty-two other trace element concentrations in the natural river water certified reference material SLRS -6 (NRC -CNRC)[END_REF] (Table S3).

The REEs concentrations were normalized by the Post Archean Australian Shale (PAAS) [START_REF] Mclennan | Rare earth elements in sedimentary rocks; influence of provenance and sedimentary processes[END_REF][START_REF] Rétif | From geochemistry to ecotoxicology of rare earth elements in aquatic environments: Diversity and uses of normalization reference materials and anomaly calculation methods[END_REF]. The Ce and Gd anomalies were calculated using Eq (1)

and Eq(2) [START_REF] Hissler | Anthropogenic rare earth element fluxes into floodplains: Coupling between geochemical monitoring and hydrodynamic sediment transport modelling[END_REF], respectively.

Ce/Ce*=CePAAS/(0.5×LaPAAS+0.5×PrPAAS) (1)

Gd/Gd* = GdPAAS/(0.4xNdPAAS + 0.6xDyPAAS) (2)
where X is the measured concentration in samples, X* the calculated geogenic concentration, and XPAAS the concentration normalized by the PAAS. Finally the anthropogenic Gd concentration was deduced using Eq(3):

Gdanth = Gd -Gd* (3)
It should be noticed that the calculation method used in the present study for the Gd anomaly tends to maximize the value of the Gd anomaly compared with other calculation methods (Louis et al., 2020;[START_REF] Rétif | From geochemistry to ecotoxicology of rare earth elements in aquatic environments: Diversity and uses of normalization reference materials and anomaly calculation methods[END_REF].

Optical methods

Synchronous fluorescence spectra (SF50) were collected on a Hitachi F-2500 fluorimeter equipped with a Xenon lamp, by using FL Solution 2.0 software and a 1 cm x 1 cm 3.5 mL quartz cuvette. The gap between excitation and emission was set at 50 nm. The 50 nm deviation provides information on both fluorescence due to protein-like substances (excitation at around 280 nm) [START_REF] Baker | Fluorescence properties of some farm wastes: implications for water quality monitoring[END_REF] and humic substances (humic acids and fulvic acids;

excitation between 300 nm and 400 nm) [START_REF] Coble | Characterization of marine and terrestrial DOM in seawater using excitation-emission spectroscopy[END_REF]. Since the absorbance of the sample may interfere with fluorescence, the absorbance at 254 nm was checked for each sample. UVvis spectra (200-600 nm) were collected on a Shimadzu UV-2600 spectrophotometer, using a 1 cm x 1 cm quartz 3.5 mL cuvette. Ultra-pure water was used for blanks.

As A254 was less than 0.1 cm -1 for 90% of the samples, the absorbance was considered low enough to record fluorescence directly, without any correction for the inner-filter effect. The SF50 spectra were collected at the natural pH of the samples. The blank was performed with ultra-pure water. The fluorescence intensities were expressed in Raman units (R.u.).

To extract quantitative information out of the SF50 spectra, a decomposition procedure was applied [START_REF] Assaad | Spectrophotometric characterization of dissolved organic matter in a rural watershed: the Madon River (N-E France)[END_REF]. In this approach, the synchronous fluorescence spectrum of a single fluorophore is represented by a Gauss function (Eq (4)):

𝑆𝐹() = 𝐹(𝑏) • 𝑒𝑥𝑝(-( -𝑏) 2 2𝑐 2 ⁄ ) (4

Databases

Maps were constructed using Qgis 3.16 (https://www.qgis.org). Shape files for the Danube River and its main tributaries were extracted from the WISE database (https://www.eea.europa.eu/data-and-maps/data/wise-large-rivers-and-large-lakes).

Brigach and Breg courses were obtained from the Baden-Württemberg Government (https://rips-dienste.lubw.baden-wuerttemberg.de).

The location and the capacity of the urban wastewater treatment plants along the Danube River and its (sub)tributaries were obtained from the dissemination platform of the Urban Waste Water Treatment Directive (UWWTD) (https://uwwtd.eu/) for the EU countries (2016 data) and from ICPDR for non-EU countries.

The 1:1,000,000 geological (GK1000) data of Germany were downloaded from the database of the federal Institute for Geosciences and Natural Resources (https://www.bgr.bund.de).

Ancillary data for DOC and nitrates were compared to the 2017 water quality data provided

Results and discussion

Physical-chemical parameters

This section summarizes the main characteristics of the surface water of the Danube, and particularly those related to dissolved organic matter (DOM) and nitrates which can provide information on the global anthropogenic pressure on the river basin.

In their upper parts, the Brigach River and the Breg River are flowing on silica-rich rocks (granite, sandstone). Their concentrations in alkaline minerals (calcium, magnesium, potassium and sodium) increase gradually moving downstream where they are leaching Triassic and Jurassic rocks (i.e. limestone, marlstone) (Figure S1 and Table S5). At Ulm, the Danube River meets the Iller River, its first large tributary (Figure 3). This alpine river average flowrate at the confluence is larger than the Danube River average flowrate (70 m 3 /s versus 53 m 3 /s), which causes a strong dilution. This dilution (Figure S2) is further accentuated at the junction with the Inn River, a large alpine tributary of low mineral content (Table S5) whose discharge flowrate is similar to that of the Danube at their junction (≈ 2200 rkm). The Inn River has its maximal flowrate in summer due to the snow melting in the Alpine upper part of its watershed (Figure S3). Its dilution effect is seen on the mineral content of the Danube waters (Figure S2). Downstream Passau, the calcium, magnesium and potassium concentrations are stable until the Danube reaches the Black Sea. The sodium concentration increases, which can be seen as a result of the discharge, in the aquatic environment, of urban wastewater containing detergent residues. DOC is higher in the samples collected upstream Ulm, in particular in the waters of the Brigach River and of the Danube River, than downstream Ulm (average of 3.1 mg /L and 2.0 mg /L respectively) (Figure 4). The Inn River carries very little DOC (0.5 mg/L in Passau-Ingling) (Table S4). DOC does not vary much in the middle section of the Danube River (average of 1.7 mg/L).

Large inflows from the tributaries (up to 7.5 mg /L for Morava River) are quickly diluted: the Morava River discharge rate is only 6% of the Danube River at their junction. In the Lower Danube, DOC increases (average 2.8 mg /L), which seems to be controlled by inflows from tributaries with less dilution effects than upstream. This increase in DOC in the Lower Danube was also highlighted during the JDS3 campaign (Joint Danube Survey 3, 2015). Synchronous fluorescence spectra and their deconvolution provide further insight into the nature of DOM (Figure 5). Protein-type fluorophores are linked to biological reactions in water but also to the runoff of biological substances from the watershed and the discharge of raw sewage or insufficiently treated urban wastewater. Their fluorescence is largely due to the presence of the indole group which is present in tryptophan (an essential amino acid for humans, which is produced only by microorganisms and plants), in auxins (plant hormones)

and in acids bound to urine (such as, for example, 5-hydroxy-indolylacetic acid, a metabolite of serotonin). Other fluorescent substances that can be found in the same spectral region as tryptophan in the SF50 spectra are threonine (λexc ≈ 280 nm), another essential amino acid for humans, and tyrosine (λexc ≈ 270 nm), also an amino acid.

In the Danube watershed, the total pseudo-concentrations of fluorophores vary from 0.042 to 0.235 Raman units. These pseudo-concentrations are globally higher for the samples from the Brigach River and the Upper Danube River (0.065 to 0.235 Raman units). They decrease in the Median Basin (0.042 to 0.105 Raman units) before increasing after the confluence with the Timok (SW42-Tim) in the Lower Basin and the Delta (0.074 to 0.141 Raman units). The Eisenbach River, the Iller River and the Inn River carry little fluorescence. Other tributaries such as the Morava River, the Hron River, the Ipel River, the Tisza River and the Sava River bring a high quantity of fluorescent dissolved organic matter in the Upper and Middle sections of the Danube River Basin (0.114 to 0.208 Raman units). However, these contributions seem to disappear rapidly under the effects of dilution with the Danube River waters, as observed for the DOC values. In the Lower section of the Danube River basin and in its delta, the contribution of other tributaries (Iskar River, Yantra River, Russenski Lom River and Prut River) is also significant (0.103 to 0.211 Raman units) and seems less affected by dilution effects, as the fluorescence of the River Danube waters is higher downstream. In conclusion, once the Danube River meets its first large tributaries (Iller River in Ulm and Inn

River in Passau) the potential anthropogenic pressure exerted by the watershed is rather diffuse and no hot spot can be detected, either for nitrogen species or dissolved organic carbon. In particular, the dissolved organic carbon is slowly increasing, which can be explained by the increase of the served population from upstream to downstream.

REE spectra and anomalies

Figure 7 shows the REE normalized patterns for the Breg River and the Brigach River. The most upstream section of both rivers present patterns with flat sections between Dy and Lu (Heavy REE or HREE). These patterns are similar to those observed by Louis et al. (2020) for headwaters in the Vosges mountains, running on a similar geology (granite, sandstone). A negative Ce anomaly, induced by the redox characteristics of cerium, can be noticed for both rivers [START_REF] Elderfield | The rare earth elements in rivers, estuaries, and coastal seas and their significance to the composition of ocean waters[END_REF]. Indeed, Ce can have two oxidation states (Ce 3+ and Ce 4+ ) in the typical conditions of temperature and pressure of surface water and Ce 4+ is likely to combine with oxygen to form cerianite. As cerianite is very insoluble, its formation leads to Ce negative anomaly in water [START_REF] Seto | Chemical condition for the appearance of a negative Ce anomaly in stream waters and groundwaters[END_REF]. Along the Danube River, all samples but two exhibit a negative Ce anomaly (Figure S5) with an average value (± one standard deviation) of 0.66 ± 0.14 between Donaueschingen and Ulm and a slight trend downstream of Ulm (coefficient of determination R 2 = 0.14).

The La/Lu ratio varies between 0.05 and 0.35 in the Breg and Brigach river and tends to increase along their courses indicating a gradual depletion in LREE with respect to HREE (Figure S6). No clear trend is observed between Donaueshingen and Ulm for the La/Lu ratio because of the variable geology (marlstone, limestone) of the Danube basin in this section combined with possible anthropogenic influences (see below). The sampling is not sufficiently dense in this section to go deeper in the discussion. Between Passau (junction with the Inn River) and the Black Sea, the La/Lu ratio decreases slowly (R 2 = 0.18), revealing an increasing fractionation of the REE pattern. The observed decrease indicates increasing inputs from urban and industrial wastewater treatment plants along the Danube River, in agreement with the increase in sodium concentration. Indeed, highly fractionated REE patterns with LREE depletion and marked Gd anomalies are typical of wastewaters [START_REF] Martin | Contrasting distribution of REE and yttrium among particulate, colloidal and dissolved fractions during low and high flows in peri-urban and agricultural river systems[END_REF]. The explored tributaries present La/Lu ratios lower than those observed on the Danube River, suggesting a more marked impact of wastewater discharges in these water courses.

Gd anomalies are detected on the Breg River at the Vörenbach station (Be4) and on the Brigach River at the Sankt Georgen station (Bi3), about 10 km downstream their sources. All the samples collected between Donaueschingen and Ulm present a strong Gd anomaly (between 2.8 at station D2 and 10.9 at station D12). S4). To define a Gd anomaly as an anthropogenic anomaly, the threshold value of Gd/Gd* is usually set to 1.5 (Bau et al., 2006). In the case of the Danube River and its tributaries, the visible Gd anomalies are always above this threshold value: that allows the calculation of Gdanth concentrations. For the Danube River itself, the most important Gd/Gd* value is for the sample JDS4-4 (Gd/Gd*: 25.9). This sampling station is located in Germany, just downstream a WWTP effluent discharge.

Only one other study has reported a positive Gd anomaly in Danube River. This is by Kulaksız and Bau (2011b), who analysed a sample grabbed in Austria. Using their data with the Nd and Dy calculation method, the Gd/Gd* anomaly is 2.8 and the Gdanth concentration is estimated at 47 ng/L. The JDS4 samples from SW6-D to SW10-D were collected in Austria. Their Gd/Gd* values are higher (from 10.8 to 18.1) for a lower calculated Gdanth concentration (from 8.3 to 9.5 ng Gdanth/L), which can be explained by the fact that the overall REE amounts for JDS4 are lower than for the study by Kulaksız and Bau, (2011b). However, as the exact location and sampling period of the Kulaksız and Bau (2011b) study are not known and it is difficult to establish a meaningful comparison.

Based on the values of Gd/Gd* and the corresponding calculated Gdanth concentrations two groups can be identified clearly along the Danube River and in its tributaries: from sample SW-D1 (Boefinger Halede) to sample SW-D31 (Ilok), with a Gd/Gd* globally above 10 and from sample SW-32 (Tisza) to sample SW-51D (Vilkova), with a Gd/Gd* globally below 10. Gd anomaly and Gdanth concentration in these two groups can be related to the number of MRI units and MRI exams in each crossed country (Fig. 10). It is interesting to note that the average Gdanth concentration between SW1 and SW31 (11 ng/L ± 2.5 n/L) is quite stable and does not closely follow the fluctuations in Gd anomalies. In the first group (Austria, Czech Republic, Slovakia, Hungary, Croatia), the number of MRI exams/year/inhabitants is between 4,000 and 6,000, while it exceeds 14,000 exams/year/inhabitants in Germany (Figure S7). In the second group (Serbia, Bulgaria, Romania), MRI exams/year/inhabitants are below 2,000 and data are missing for Moldavia and Ukraine. The low number of MRI exams in Bulgaria can explain why there is no Gd/Gd* detected in sample JDS4-46. The tipping point, where Gd/Gd* becomes lower than 10, is Ilok (Croatia) (i.e. SW31-D). We note that the sampling plan of JDS4 for the Danube Delta region included only one site in the Chilia arm. The possible impact of the large city of Tulcea on REE pattern in the Sulina and Sfantu Gheorghe arms of the Danube Delta would require further consideration.

The groundwater wells close to the Danube River are also affected by anthropogenic Gd (Fig. S8). Gd-containing contrasting agents are synthetized not to easily release ionic Gd 3+ , which is extremely toxic, inside the human body. This otherwise desirable property makes them particularly difficult to degrade during drinking water treatment production and explains the presence of Gd in tap water.

It can be expected that the situation will not improve in the near future [START_REF] Lachaux | Integrated environmental risk assessment of rare earth elements mixture on aquatic ecosystems[END_REF].

The risk associated with Gd might even increase with the application of advanced technologies The annual flux of Gdanth to the Black Sea, based on a total discharge just upstream of the delta of 6,500 m 3 /sec, is estimated at 0.9 ton/year, corresponding to 11 kg/year/million inhabitants. This is of the same order of magnitude as the flux discharged by the Garonne River in France (14 kg/year/million inhabitants) (Lerat-Hardy et al., 2019). The annual flux is comparable with the annual flux of other large watersheds, such as the Saint-Laurent Estuary and Gulf in Canada (0.9 ± 0.7 ton/year) [START_REF] Dang | The contrasting estuarine geochemistry of rare earth elements between ice-covered and ice-free conditions[END_REF] and the Pearl River (0.6 ± 0.6 ton/year for about 84 million inhabitants) (Ma and Wang, 2023).

Finally, an Ytterbium anomaly seems to appear at some sampling sites close to the delta, in particular in stations SW-D37 and SW-D48. This anomaly is also clear on two tributaries, the Sava River and the Velika Morava River. No information related to such an anomaly could be found in the literature. Ytterbium is not largely used: its main uses are related to the improvement of stainless-steel properties and the production of lasers. However, sedimentary phosphate deposits (Emsbo et al., 2015) and their dumped residues from fertilizer production can contain ytterbium [START_REF] Zirnea | Cluster analysis based of geochemical properties of phosphogypsum dump located near Bacau City in Romania[END_REF], which can be leached toward water bodies.

Conclusions

At between Mohringen (Germany) and Ulm. Codes in Table SI4.

Figure 1 :

 1 Figure 1: Sampling sites along the Danube River between Ulm and Donaueschingen (samples

Figure 2 :

 2 Figure 2: Sampling sites along the Danube River between Ulm and the Black Sea (green dots)

  )where F(b) is the pseudo-concentration of the fluorophore, b (in nm) its characteristic excitation wavelength, i.e. the wavelength where the maximal fluorescence is recorded and c a parameter related to the width of the Gauss function. Each SF50 spectrum was decomposed into five Gauss functions (i.e., distinct fluorophores groups) by deconvolution. With a 50 nm gap the water Raman scatter interferes with the fluorescence of the fluorophore with a characteristic excitation wavelength of 355 nm and the contribution of water to this fluorophore is further corrected for.

Figure 3 :

 3 Figure 3: Danube River watershed area (in orange) as a function of the received tributaries

Figure 4 :

 4 Figure 4: Dissolved organic carbon (DOC) along the Danube River. D = Donaueschingen, U =

Figure 5 :

 5 Figure 5: Protein-like and humic-like fluorescence in the Danube River watershed in July 2019.

Figure 6 :

 6 Figure 6: Nitrates along the Danube River. D = Donaueschingen, U = Ulm.

Figure 7 :

 7 Figure 7: REE patterns, normalized with respect to PAAS, along the Breg River and the Brigach

Figure 9 :

 9 Figure 9: Gd anomaly and Gdanth concentration downstream from the sources to Ulm (a) and

Figure 10

 10 Figure10was built using data on the location and capacity of wastewater treatments plants

  designed to degrade urban micropollutants in wastewater treatments plants, with the ultimate goals to reduce aquatic ecosystems contamination on one hand and potentially reuse treated wastewater on another hand. The risk associated with the dissemination of free Gd 3+ cannot be ignored considering the possible degradation of Gd-based contrast agents(Oluwasala et al., 2022).

Figure 10 :

 10 Figure 10: Gd anomaly in the Danube River watershed in as a function of the population and

  short distance from their sources the Brigach River and Breg River present REE patterns already influenced by the release in the aquatic environment of anthropogenic Gd. The Gd concentration in the Danube River is influenced by the development of MRI facilities and examinations. The slow decrease of the La/Lu ratio between Ulm and the Black Sea indicates an increasing anthropogenic influence as the Danube River flows downstream. As this survey was the first one to analyze REEs in the Danube River basin, future work will focus on consolidating the data with an increase of the number of sampling points in the Middle and Lower sections of the Danube as well as along the largest tributaries. Analysis of sediments could also bring information on the distribution of REE in the watershed. Although the monitoring of REE and the understanding of their origin and fate in large and complex watersheds are not easy tasks, they should be pursued, taking the example of the logistics developed during JDS4.

Figure S1 :

 S1 Figure S1: Main geological characteristics a) between the sources and Nendingen (Germany) b)

Figure

  Figure S2: Main anions concentrations along the Danube River during JDS4. Comparison with the TNMN data for 2017. Donaueschingen = D, Ulm = U

  Flow rates in some specific stations were measured during the JDS4 and retrieved from the JDS4 database (https://data.danubesurvey.org/jds4/).The Eurostat database provided the data about the number of MRI examinations in Europe (https://data.europa.eu). Complementary data for non-EU countries were extracted from the OECD Health Database (https://www.oecd.org/els/health-systems/health-data.htm).The distances are established from the mouth at the Black Sea for the Danube River and from the confluence with the Danube River for its tributaries. The Danube River is divided into three main parts: the Upper Danube River refers to the river course between the sources and rkm 1789.5 (Kližská Nemá, Slovakia); the Middle Danube River to the river course between rkm 1789.5 and rkm 943 (Turnu Severin, Romania); the Lower Danube River to the river course

	between	rkm	943	and	the	Black	Sea

by ICPDR through the TransNational Monitoring Network (TNMN) scheme (https://www.icpdr.org/main/sites/default/files/nodes/documents/tnmn2017_annex_1.xls).

(https://www.icpdr.org/flowpaper/viewer/default/files/roof_report_2004_annexes.pdf).

Acknowledgments

The authors acknowledge the financial support from the National Agency for Research -France (ANR) through the project ANR-16-CE34-0012-001 and of the LIFE Program (ENV.780956). The support of ICPDR and Norman network is greatly acknowledged.

Spatial distribution of Rare Earth Elements in a transnational watershed: the case of the Danube River

Pauline LOUIS 1 , Davide A.L. VIGNATI 2 , Steve PONTVIANNE 1 , Marie-Noëlle PONS 1,3 * 1 Université de Lorraine, CNRS, LRGP, F-54000 Nancy, France 2 Université de Lorraine, CNRS, LIEC, F-57000 Metz, France 3 LTSER-Zone Atelier du Bassin de la Moselle, LRGP, 54000, Nancy, France Supplementary Information


MD Grégory Verhoest
email: gregory.verhoest@wanadoo.fr

Raffi Avakian

Karim Bensalah

Rodolphe Thuret

Vincenzo Ficarra

Walter Artibani

Jacques Tostain

François Guille

Lucas Cindolo

Grégory Verhoest¹

Raffi Avakian¹

Karim Bensalah¹

Vincenzo Ficarra³

Walter Artibani³

François Guille¹

Alexandre De

L A Taille

Claude-Clément Abbou

Laurent Salomon

Nathalie Rioux-Leclercq

Jean-Jacques Patard¹

Urinary collecting system invasion is an independent prognostic factor of organ confined renal

Keywords: kidney cancer, collecting system, tumor invasion, pathology, nephrectomy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction:

Renal cell carcinoma (RCC) represents 2-3% of all solid neoplasms and is the third most common genito-urinary tumor. Clear cell histological subtype accounts for 85% of the RCC cases [START_REF] Pantuck | The changing natural history of renal cell carcinoma[END_REF] while the annual incidence of renal tumors is increasing in by around 3% (1.5-5%). This may be due in part to the widespread use of imaging techniques such as ultrasound (US) and computerized tomography (CT) [START_REF] Patard | The changing evolution of renal tumours: a single center experience over a two-decade period[END_REF] which allows the early detection of an increasing number of small renal masses.

The value of an accurate cancer staging system is crucial, not only for defining specific prognostic tools that will identify groups with greater risk of recurrence but also for providing cost-effective post-therapeutic surveillance protocols. TNM staging system remains the most useful determinant in disease classification reflecting the propensity of these tumors to extend locally (T), involve regional lymph nodes (N) and initiate metastatic growth (M) [START_REF] Ficarra | Multiinstitutional European validation of the 2002 TNM staging system in conventional and papillary localized renal cell carcinoma[END_REF] .

Although invasion of Gerota's fascia, perirenal fat, adrenal gland, renal vein, inferior vena cava, and adjacent retroperitoneal structures have been extensively studied and discussed as prognostic variables, the prognostic value of the urinary collecting system invasion (UCSI) has not focused great interest in the recent RCC prognostication literature [START_REF] Ficarra | Proposal for revision of the TNM classification system for renal cell carcinoma[END_REF] . Invasion of the urinary collecting system, mentioned in the 1978 edition of the TNM classification, is no longer considered as a criterion for staging RCC [START_REF] Sobin | TNM. Classification of malignant tumors[END_REF] . However, some recent studies suggested that UCSI could be potentially interesting as a prognostic variable in low stage tumors [START_REF] Ficarra | TNM staging system for renal-cell carcinoma: current status and future perspectives[END_REF] .

The aim of this large multicenter retrospective study was to assess the relationship between UCSI and usual prognostic variables in RCC and to evaluate to what extent UCSI should be considered as a major prognostic variable in RCC leading to further change of the TNM classification.

Material and Methods

Data retrieval

This retrospective study included 1124 patients who underwent either radical nephrectomy or nephron sparing surgery for a renal tumor between 1997 and 2004. Data were extracted from RCC databases of 5 European academic centres: Rennes, Créteil, Saint-Etienne (France), Padua and Benevento (Italy). The study was approved by the IRB at each institution. Only tumors with malignant histology at final pathology were retained for analysis. Patients with metastatic disease were treated by radical nephrectomy followed by interferon, interleukin 2based immunotherapy according to the standards of care at that time. Post-operative follow-up was established according to protocols at each institution.

Clinical and histological variables.

Involvement of the urinary collecting system (UCS) was defined as evidence of invasion of the calices, infundibulum, pelvis or ureter; excluding simple compression or displacement by the primary renal mass. Tumor stage and grade were determined according to the 2002 TNM classification [START_REF] Sobin | TNM. Classification of malignant tumors[END_REF] , and Fuhrman grading system respectively [START_REF] Fuhrman | Prognostic significance of morphologic parameters in renal cell carcinoma[END_REF] . Histologic subtypes were stratified according to the World Health Organization classification [START_REF] Lopez-Beltran | WHO classification of the renal tumors of the adults[END_REF] . All the following variables were recorded and analysed: age, gender, presence of symptoms at diagnosis, Eastern Cooperative Oncology Group (ECOG) Performance Status (PS), tumor size, TNM stage, Fuhrman grade, UCSI, histologic subtype, and survival data. At last follow up, patients were considered as living (with or without disease) or deceased (from renal cancer or from any other cause).

Statistical analysis:

For statistical analysis, patients were classified as ECOG 0 versus 1 or higher. Qualitative and quantitative variables were compared by using Chi-square (Fischer exact test) and Student t tests, respectively. All analyses were conducted with the statistical Package for the Social Sciences version 13.1 (SPSS Inc, Chicago, IL, USA), and p value significance was fixed at 0.05. Univariate and multivariate analyses were performed by using the Kaplan-Meier and the Cox methods respectively.

Results

Patients and tumor characteristics

This study included 1124 patients; out of them 771 were males (68.6%) and 353 were females (31.4%). Median age was 61 years [14-88]. Symptoms were present at diagnosis in 50.3% of the cases. Median tumor size was 6 cm [1-24]. Tumors were organ confined and recorded as Fuhrman grade 1 or 2 in 67.2% and 62.3% of the cases respectively. Surgical procedure was a radical nephrectomy in 896 cases (79.7%). The most common histopathological diagnosis was clear cell carcinoma in 875 cases (84%). Mean follow-up was 63 months. At the end of the follow-up, 246 patients (21.9%) died from cancer. Patient and tumor characteristics are summarized in Table 1. One hundred and thirty two patients (11.7%) had pathological evidence of UCSI.

Relationship between UCSI and clinico-pathological features

UCSI was clearly associated with aggressive clinico-pathological features. For example, when symptoms were present at diagnosis, UCSI was observed in 17.6% of the cases vs only 4.6% in incidental tumors. Similarly, when ECOC PS was recorded as 1 or more, 17.4% of the tumors invaded UCS vs 9.5% in case of ECOG 0. Locally advanced (T Stage ≥3) and high grade tumors (Fuhrman grade ≥3) were associated with more frequent UCSI (26.2% and 23.1 %) compared to organ confined and low grade tumors (4.6% in both cases, p <0.001). In case of nodal invasion or distant metastases, UCSI was also more frequently noted (34.3% and 22.3%, p <0.001). However, histological subtype had no impact on UCSI. All results regarding relationship between UCSI and clinico-pathological features are summarized in table 2.

Survival analysis

Evaluation of survival outcome of patients without UCSI revealed that cancer specific survival at 5 and 10 years was 80% and 74% respectively, while patients with UCSI had 43% and 41% cancer specific survival rates at 5 and 10 years respectively (Figure 1). Median survival time was not reached for patients without UCSI, whereas it was 60 months when the collecting system was invaded (p = 0.0001). Sub stratifying T stage with or without UCSI revealed a significant survival difference either for patients with pT1-T2 tumors (p = 0.0001), or for those with pT3-T4 tumors (p = 0.04) (Figure 2). When comparing pT2 tumors with UCSI to pT3 tumors there was no survival difference (p = 0.95).

In univariate analysis, TNM stage, Fuhrman grade, symptoms at diagnosis, ECOG-PS, median tumor size and UCSI (p <0.001) appeared as significant prognostic variables for cancer specific survival (Table 3).

In multivariate analysis, TNM stage, median tumor size (p <0.001), ECOG-PS (p = 0.002) and Fuhrman grade (p = 0.02) were independent prognostic factors. Interestingly, UCSI was not retained as an independent prognostic factor for the entire population (Table 3). However, when searching separately for independent prognostic factors in pT1-T2 and pT3-T4 tumors, it appeared that UCSI was an independent prognostic variable in pT1-T2 tumors while it was not in pT3-T4 tumors (tables 4 and 5). When both ECOG performance status and symptoms at diagnosis were included within the organ confined Cox model, only ECOG was retained as an independent prognostic parameter. When symptoms at diagnosis were included within the model but not ECOG performance status; both UCSI (p = 0.02) and symptoms at presentation (p = 0.002) were retained as independent prognostic parameters. Consistently with the current 2002 TNM classification, perinephritic fat invasion (p <0.001) and renal vein invasion (p = 0.01) were retained as independent prognostic factors in the setting of pT3-T4 tumors. Median tumor size was also found to be an additional independent prognostic parameter in pT3-T4 tumors (Table 6).

Discussion

There is a consensus that TNM classification needs to be continuously revised according to the current state of the art in prognostication [START_REF] Gospodarowicz | The process for continuous improvement of the TNM classification[END_REF] . Large retrospective databases allow such a comprehensive analysis of anatomic and pathologic parameters that are the basis of the TNM classification. For example, the optimal tumor size cut-off for delineating T1 from T2 tumors has been extensively discussed [START_REF] Karakiewicz | Tumor Size Improves the Accuracy of TNM Predictions in Patients with Renal Cancer[END_REF] . Similarly, T3 sub classification of the 2002 TNM classification has been extensively challenged with discussion about prognostic value of perinephritic fat invasion versus sinus fat invasion [START_REF] Thompson | IS RENAL SINUS FAT INVASION THE SAME AS PERINEPHRIC FAT INVASION FOR pT3a RENAL CELL CARCINOMA?[END_REF] , adrenal gland invasion [START_REF] Thompson | Should direct ipsilateral adrenal invasion from renal cell carcinoma be classified as pT3a?[END_REF] , renal vein and the level of IVC invasion [START_REF] Ficarra | Proposal for reclassification of the TNM staging system in patients with locally advanced (pT3-4) renal cell carcinoma according to the cancer-related outcome[END_REF] . By contrast, prognostic value of UCSI has been poorly investigated so far.

To our knowledge, our series is the largest series published in the literature, with a cohort of 1124 patients coming from 5 European academic centers and including 132 tumors invading the collecting system (11.7%). UCSI frequency that we reported is consistent with the range of the literature going from 7.5 to 14% [START_REF] Uzzo | Renal cell carcinoma invading the urinary collecting system: implications for staging[END_REF][START_REF] Palapattu | Collecting system invasion in renal cell carcinoma: impact on prognosis and future staging strategies[END_REF][START_REF] Terrone | Prognostic value of the involvement of the urinary collecting system in renal cell carcinoma[END_REF][START_REF] Klatte | Prognostic relevance of capsular involvement and collecting system invasion in stage I and II renal cell carcinoma[END_REF] . We thus demonstrated that tumors with UCSI presented with higher stages and grades and were more likely to invade lymph nodes or exhibit distant metastases. Consistently, these tumors which were larger in size were more likely to be symptomatic. Interestingly, although UCSI was not retained as an independent prognostic parameter overall, it was retained as an independent prognostic variable in pT1-T2 stage tumors. It suggests that pathological parameters such as renal vein invasion or perinephritic fat invasion are more important outcome predictors than UCSI in locally advanced tumors. Conversely, it suggests that UCSI, together with tumor size, ECOG performance status and symptoms could be an additional important prognostic variable in organ confined renal tumors. Interestingly we previously demonstrated that symptoms as well as tumor size were able to stratify T1-T2 tumors [START_REF] Patard | Symptoms as well as tumor size provide prognostic information on patients with localized renal tumors[END_REF] . This additional data could provide some rationale for using tumor size, symptoms and UCSI for sub stratifying T1-T2 tumors similarly to that of fat, adrenal and renal vein invasion which are defining T3a, b and c stages in locally advanced tumors. We and others also demonstrated that tumors associated with symptoms at diagnosis were larger in size, of higher grades and exhibited poorer outcome [START_REF] Patard | Multi-institutional validation of a symptom based classification for renal cell carcinoma[END_REF] . Therefore, local symptoms such as haematuria could be viewed as surrogate indicators for UCSI.

However, interestingly the present study highlighted the independent prognostic value of pathological invasion of UCS even after adjustment for tumor size, ECOG performance status or symptoms. It means that UCSI has an intrinsic biological significance in organ confined tumors, similarly to that of renal vein or capsular invasion in locally advanced disease. A tumor crossing an anatomic natural barrier such as collecting system should not be longer considered as an organ confined tumor. This is supported by the fact that we found no survival difference between pT2 tumors exhibiting UCSI and pT3 tumors. It should therefore be taken into account for further TNM classification.

Only 4 studies have focused on UCSI prognostic significance so far. In 2002, Uzzo et al analyzed the outcome of 61 patients with tumors invading UCS in RCC [START_REF] Uzzo | Renal cell carcinoma invading the urinary collecting system: implications for staging[END_REF] . They concluded that these tumors although exhibiting higher stages and grades had not a significantly poorer outcome. It was noticed that even though UCSI in organ confined RCCs was uncommon, there was a trend for a worse prognosis when it occurred. However, the independent prognostic value of UCSI was not assessed through multivariate analyses in this study. By analyzing 124 tumors invading the UCS in 895 patients, Palapattu et al reported a high proportion of large (2/3 were classified as pT3) and symptomatic tumors within this population. UCSI was also significantly associated with lymph node invasion and presence of distant metastases at diagnosis [START_REF] Palapattu | Collecting system invasion in renal cell carcinoma: impact on prognosis and future staging strategies[END_REF] . Univariate analysis revealed a 1.96 time increased risk of death for patients with UCSI. Interestingly, the authors noticed a significantly worse prognosis in patients with T1 and T2 tumors with UCSI. Recently, the same group analyzed the outcome of 39 patients with UCSI in stage I and II RCCs [START_REF] Klatte | Prognostic relevance of capsular involvement and collecting system invasion in stage I and II renal cell carcinoma[END_REF] . No association was found between UCSI, T stage, Fuhrman grade, tumor size, ECOG-PS, capsular involvement and histological type. Univariate analysis confirmed a worse prognosis when UCS was invaded, resulting in a significantly poorer survival in this subgroup of patients. Furthermore consistently with our study, multivariate analysis revealed UCSI as an independent prognostic factor within organ confined tumors, with a 3.78 times increased risk of tumor recurrence.

Finally, Terrone et al obtained comparable results except that adverse outcome associated with UCSI was only observed within the pT2 category [START_REF] Terrone | Prognostic value of the involvement of the urinary collecting system in renal cell carcinoma[END_REF] . The current series dealing with UCSI prognostic value are summarized in table 7.

Even though TNM staging remains the most widely recognized prognostic system it is important to further identify independent prognostic variables that are able either to improve the system itself or that can be integrated within more sophisticated prognostic tools. The common goal of identifying new prognostic variables is to improve the predictive accuracy of current prognostic systems that we dispose. Based on our results, we believe that UCSI could fit with such an objective for organ confined tumors. Based on the current TNM classification, organ confined renal tumors are only stratified based on tumor size. No reliable biological prognostic variable has been indentified so far in the setting of localized RCC. Therefore, information regarding UCSI could be either integrated within further TNM classification revision or being tested for incorporated within new nomograms specially designed for localized RCC. Finally, our study besides new insights in localized RCC prognostication also suggests that a closer follow up should be proposed for those patients with organ confined tumors a nd UCSI by contrast with what is usually proposed in such low stage RCCs 20 .

Our study is not exempt of limitations. The main drawback is that due to its retrospective character and to the lack of central pathological review, some critical information may have been missed. Notably, it could be argued that since UCSI has not been prospectively searched, only massive UCSI associated with large tumors has been reported.

Then less obvious UCSI either microscopic or limited may have been missed even though it was present, therefore leading to a selection bias of more aggressive tumors. Additionally, because our study was based on retrospective pathological reports analysis, it could be hypothesized that this information was more likely to be noted in larger tumors that had been treated by radical nephrectomy rather than in small tumours treated by nephron sparing surgery. Our low NSS rate pleads for this hypothesis. Therefore by searching for CSI information based on the current pathological analysis practice we biased for larger tumors which are more likely to be symptomatic. Our 50% symptomatic tumors rate which is a little bit higher than usually reported is also consistent with this hypothesis. Finally, it could also be argued that less attention is paid for CSI in small renal tumor specimens and that CSI is more difficult to assess following partial nephrectomy. It is probably true, however even though the size of tumor margins does not matter, careful examination of the quality of the tumor surrounding tissue is usually performed. Therefore, the presence or absence of urothelium tissue along with its invasion should be easily identified. Urologists should also mention whether or not collecting system has been entered and if there is any CSI suspicion.

Conclusion

Tumors with UCSI are associated with most adverse clinical and pathological parameters in RCC. Based on the present study, we recommend that urinary collecting system invasion should be systematically identified and recorded in pathological reports. UCSI appears as a new important prognostic factor for organ confined tumors. This pathologic parameter should be considered for further TNM classification revisions and for its integration within new prognostic nomograms. Finally, a closer follow-up should be proposed to patients operated on for organ confined tumors and exhibiting UCSI. To what extend organ confined tumors with associated UCSI should be considered as locally advanced tumors and proposed for adjuvant trials needs to be further explored.