R Ceolin

S Toscani

Ivo B Rietveld

M Barrio

J Ll

Tamarit

Pitfalls and feedback when constructing topological pressuretemperature phase diagrams

The stability hierarchy between different phases of a chemical compound can be accurately reproduced in a topological phase diagram. This type of phase diagrams may appear to be the result of simple extrapolations, however, experimental complications quickly increase in the case of crystalline trimorphism (and higher order polymorphism). To ensure the accurate positioning of stable phase domains, a topological phase diagram needs to be consistent. This paper gives an example of how thermodynamic feedback can be used in the topological construction of phase diagrams to ensure overall consistency in a phase diagram based on the case of piracetam crystalline trimorphism.

Introduction

In a recent article, a pressure-temperature phase diagrams for the trimorphism of piracetam has been obtained through a topological approach [1]. An annotated version of the stable pressure-temperature (P-T) phase diagram proposed in Figure 4 of reference [1] is provided here in Figure 1. It appears that the final phase diagram of the quoted article is not in accordance with the data reported therein. For topological constructions of phase diagrams, feedback systems can be used, which can help identify and avoid inconsistencies, even if the process of constructing phase diagrams can be complicated. In the present paper, an example will be given of how inconsistencies in a topological phase diagram can be identified. 4 of reference [1]. The stable phase regions for forms I, II, and III, are indicated in accordance with ref. [1] as are the two-phase equilibria with line a: I-II, line b: I-III, and line c: II-III and the resulting triple points.

The construction of the topological phase diagram depicted in Figure 1 was based on experimental differential scanning calorimetry (DSC) and X-ray diffraction data reported in the Tables 3,4, and 7 of reference [1]. In addition, the temperature and enthalpy change of the transition II→liq was provided in the text in reference [1]. These data have been summarized in Table 1 in the present paper. For the solid-solid transition between forms III and II, neither the transition enthalpy nor the transition temperature were reported in the original publication; however, according to Figure 1, T III→II should be higher than T III→I by simple extrapolation of the two-phase equilibrium II-III to the temperature axis (P = 0). The II-III equilibrium has a steeper slope in Figure 1 than equilibrium III-I and the two equilibria cross at increased positive pressure, thus by deduction from Figure 1 the following transition temperature sequence can be obtained: T II→I < T III→I < T III→II .

The diagram in Figure 4 of reference [1] (here Figure 1) has been constructed using the alternation rule and the Clapeyron equation: (1) with dP/dT the slope of a two-phase equilibrium in the pressure-temperature phase diagram, ΔS the entropy change associated with the equilibrium, ΔV the volume change, ΔH the enthalpy change, and Τ the temperature at which the two-phase equilibrium has been observed. The alternation rule, which is a direct result of the intersecting Gibbs-energy surfaces for the different phases, states that equilibria passing through a triple point will change their stability ranking with one level and it implies that the stability ranking of the phase equilibria around a triple point will alternate between two levels of stability [2].

For the topological construction of P-T phase diagrams, data obtained at ordinary pressure (i.e. pressure conditions in a closed DSC pan) are generally used in the Clapeyron equation and equilibrium lines are approximated by straight lines. In particular solid-solid equilibria are straight over extensive temperature and pressure ranges and using this approach in the correct way leads to reliable P-T phase diagrams, as was demonstrated in the case of paracetamol [3,4], benzocaine [10], adamantane derivatives [5][6][7] and other organic compounds [8][9][10][11]. The inequalities between the slopes of the different equilibria can be determined once the respective enthalpy changes, volume changes, and equilibrium temperatures have been determined. Making use of the alternation rule, the stability ranking of the different equilibria in the phase diagram can be determined. For trimorphism, this has been demonstrated for ferrocene by two of the present authors in 2002 [12]. 4 in reference [1]. b The temperatures in kelvin have been obtained by adding 273 to the temperatures in degrees Celsius as in Table 7 of reference [1]. c By DSC as indicated in Table 4 of reference [1].

d By X-ray diffraction as indicated in Table 4 of reference [1]. e By calculation, page 170 in [1] f From Table 7 in reference [1] g Calculated by dividing the values in the fifth column by the molar weight of piracetam 142.156 g mol -1 . h From Table 7 in reference [1], and verified by applying equation 1.

On the interpretation of the calorimetric data

The topological construction of P-T phase diagrams depends for a large part on the interpretation of data obtained by DSC and X-ray diffraction under ordinary conditions. Thus, the applied data will be considered before discussing the actual construction of the phase diagram. Concerning the DSC results, it is mentioned on page 168 of the paper of reference [1]: "Randomly, for some experiments carried out by DSC on Form III at 10 °C min -1 , we get a transformation of Form III at a higher temperature, i.e. 127.5 °C (figure 1,curve b). This could not be ascribed to a transformation into Form I, the latter occurring at a lower temperature. We therefore hypothesized that this transformation could be the Form III toward Form II transformation, followed by the melting of Form II and recrystallization into Form I."

The hypothesis is incorrect for several reasons. Firstly, no enthalpy data are provided for the transition discussed in the citation above, which presumably is from form III to form II. Using a Hess cycle and starting from form III, it is a thermodynamic fact that ΔH III→I + ΔH I→II + ΔH II→III = 0 and thus that ΔH III→II = ΔH III→I + ΔH I→II (The direction of the transition II→III or III→II defines the sign of the enthalpy change). Thus with the data in Table 1, one finds ΔH III→II = 26.8 -24.8 = 2 J g -1 . It implies that the peak for the ΔH III→II must be very small in relation to those of the two other solid-solid transitions. Considering curve 'b' in Figure 1 in reference [1] (see the annotated To verify the size of the peak in form 'b', the curves have been enlarged and the sizes of the three solid-solid peaks relative to the melting peak of form I, which in each curve is the form that melts, have been determined by weighing (see Appendix). Thus, if ΔH I→liq = 187 J g -1 (taken as 100 %), a relative peak size of 15 % in curve 'a' for ΔH III→I equals 28 J g -1 and a relative peak size of 12 % in curve 'e' for ΔH II→I equals 22.4 J g -1 . The enthalpy of the unknown transition, ΔH III→x , has a relative surface of 15 % and equals therefore 28 J g -1 ; hence transition III→x is most likely III→I and certainly cannot be III→II as it was originally assigned in reference [1], because the enthalpy represented by the peak in curve 'b' is far too large for the ΔH III→II = 2 J g -1 of the III→II transition.

Secondly, once the solid-solid peak in curve 'b' in Figure 1 of reference [1] had been assigned to the III→II transition, it was assumed by the authors that the small endoexothermic effect observed a few degrees higher was the fusion of form II and recrystallization into form I, which then melts at its normal temperature. However, if form II melts and the melt turns into form I, the overall enthalpy will have to be equal to the transition enthalpy of II→I as in thermodynamics the path is not important, only the initial and the final state. Neglecting heat capacity differences, one would expect a positive enthalpy of about 24.8 J g -1 for the sum over the processes of melting and recrystallization. Determining the area of the endo-exo peak relative to the melting peak of form I, one finds a value of about 1%, which is equal to an enthalpy of 1.9 J g -1 . This cannot coincide with the proposed transitions of II→liq→I, because the enthalpy of the summed peaks is much too low.

A last point to be made is that with the data in Table 1, the slope of the III-II equilibrium can be calculated. The enthalpy change ΔH III→II was found to be 2 J g -1 . In the same way, by thermodynamic cycle, the volume change from III to II can be calculated as ∆V III→I + ∆V I→II + ∆V II→III = 0 and thus ∆V III→II ≈ 0.0103 cm 3 g -1 using the data in Table 1 (again it has to be kept in mind that the direction of the arrow defines the sign of the volume change). Using the enthalpy and volume changes together with the transition temperature of III→II of 400.5 K (127.5 °C as observed in reference [1]; see the quotation above), the slope of the equilibrium is found to be 0.48 MPa K -1 , which is a shallower slope than those of the two other equilibria. This, however, is in contradiction with the proposed phase diagram in Figure 1, where equilibrium III-II has the steepest slope. One ought to conclude therefore that the phase diagram is inconsistent with the available data. This will be discussed in more detail below.

3.-Construction of the P-T diagram using the data in Table 1

Once the data obtained by DSC and X-ray diffraction has been verified, the subsequent step is tracing the equilibrium lines using the Clapeyron equation (1). At the intersections of the equilibrium lines, the positions of the triple points can be found, which make up the framework of a topological P-T phase diagram. In the case of trimorphism, the number of triple points in a phase diagram incorporating the liquid and the vapor phase amounts to 10 as established in 1890 by Riecke [13]. For reference, the P-T phase diagram containing all 10 triple points provided in Figure 3 of reference [1], has been reproduced in Figure 2A. In Figure 2B, the phase diagram is provided that is based on the data in Table 1 and constructed in accordance with thermodynamics.

3.1.-Triple points I-II-liq, I-III-liq and II-III-liq

The three triple points I-II-liq, I-III-liq, and II-III-liq have the liquid phase in common implying that two melting equilibria and one solid-solid equilibrium intersect each other in each of the triple points. As stated in the introduction, two-phase equilibria can be safely considered straight over a considerable pressure and temperature range. Thus, using the Clapeyron equation and starting at the measured or otherwise obtained melting points in Table 1, the equations describing the melting equilibria for each polymorph can be obtained:

The pressure of each melting point (at ordinary pressure) can be safely written as 0 MPa, because the equilibrium lines are very steep relative to the pressure and the error over the values of the vapor pressures, which themselves are in the order of pascal, by taking 0 MPa will be negligible. Setting eqs. 2 and 3 equal leads to triple point I-II-liq, eqs. 2 and 4 lead to I-III-liq, and eqs. 3 and 4 lead to II-III-liq. The pressure and temperature coordinates can be found in Table 2 and their schematic positions are given in Figures 2A andB.

3.2.-Triple point I-II-III and equilibrium curve II-III obtained by triple points

The triple point I-II-III only involves solid phases. In the same way as for the solidliquid equilibria, the expressions for the equilibrium curves I-II and I-III, intersecting in I-II-III, can be obtained using the data in Table 1:

P I-II (/MPa) = 2.08 T (/K) -749 (5)
P I-III (/MPa) = 1.67 T (/K) -618 (6) Setting eqs. 5 and 6 equal leads to the (P,T) coordinates of the I-II-III triple point of 320 K and -84 MPa. The values have been compiled in Table 2 together with the data on the triple points containing the melt. Schematically, the triple point has been placed in Figure 2B as well.

It can be observed that the I-II-III triple point is located in the negative-pressure or expanded-phase domain, i.e. the triple point I-II-III must be metastable (Table 2, Figure 2B).

Considering the data in Table 1, this is not surprising as the slope of equilibrium II→I from its transition temperature of about 90°C is steeper with 2.08 MPa K -1 than the slope of equilibrium III→I with 1.67, which has a transition temperature that is located at a higher temperature of about 100 °C. This can only mean that the two equilibria diverge with increasing pressure and that the triple point involving these solid phases must be located at negative pressure (see curve a and b in Figure 2B). However, this finding is inconsistent with the phase diagram in Figure 1, where the I-II-III triple point is located at positive pressure.

The II-III equilibrium calculated above with approximate experimental values leading to a slope of 0.48 MPa K -1 , and which would confirm a negative triple point if it indeed occurs at 127°C, can also be calculated using two triple points through which the equilibrium must pass, I-II-III and II-III-liq. It results in the following expression:

P II-III (/MPa) = 2.71 T (/K) -952 (7)
The slope of this equilibrium (line c in Figure 2B) is higher than both other solid-solid transitions in accordance with the proposed phase diagram in Figure 1, however, which of the two slopes for the II-III equilibrium should be considered the right one? It can also be seen that the sequence T II→I < T III→I < T III→II as found for Figure 1 in the introduction does not hold anymore and has become T III→II < T II→I < T III→I in Figure 2B. The strength of the topological method resides in consistency in the calculated triple points and phase equilibria, which unfortunately is lacking in this case.

3.3.-Equilibrium I-III and its triple point I-III-liq

The phase diagram in Figure 2B, which has been constructed with the data in Table 1 and the equations 2 -7 derived therefrom, highlights another problem as well. It can be seen that the I-III equilibrium (line b), which is defined by eq. 6 obtained by experimental data, cannot intersect its triple point in point k in Figure 2B (point 7 in Figure 2A). This is once again an inconsistency making it impossible to construct a thermodynamically consistent phase diagram.

4.-Discussion

Comparing the topological phase diagram constructed in the present paper using the data from Table 1 in Figure 2B with the one obtained in the original paper [1] (in the present text in Figure 2A), it can be seen that they are not the same. The experimental data compiled in Table 1 do not give rise to a positive pressure for the triple point I-II-III. The inconsistency is mentioned in the original paper on page 171 [1]: "As seen on the phase diagram, since the II-I equilibrium line passes through points 5 and 6 and the III-I equilibrium line through points 4 and 7, the slope of the II-I equilibrium line is necessary lower than that of the III-I equilibrium one. According to the calculated slopes (table 7), this should have been the opposite, although very close values. But the calculation gives the slope at the origin with the approximation that the two-phase equilibrium curve is treated as a straight line. This is not completely true and may explain the disagreement between the calculation and the thermodynamic construction."

Although the remark may appear innocent enough that a slight adjustment of a small difference in the slope may explain and solve the inconsistency, it does in fact invalidate either the experimental results or thermodynamics itself. It may be possible that the experimental error allows the inversion of the slopes, but in that case the data will need to be improved before constructing the phase diagram, as the data are not conclusive. If however, the data and therefore the slopes are correct, the reasoning in the quotation above would imply that the equilibria II-I and III-I cross each other twice, once at negative and once at positive pressure. This was already excluded by Riecke in 1890, who limited the number of triple points in a phase diagram to one per possible triplet (i.e. I-II-III, I-II-liq, etc.) [13], because two-phase equilibria are monotonically rising or descending curves that can therefore only cross once, as required by thermodynamics [2]. In addition, from our experience, experimentally obtained solid-solid and even solid-liquid equilibria are in fact straight lines over considerable pressure and temperature ranges, for example piracetam [14], cysteamine hydrochloride [15], and benzocaine [16] and many others [5][6][7][8][9][10][11][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. The strength of the topological approach lies in the consistency of the results from different approaches. When there is no consistency in a topological phase diagram, the thermodynamic basis for that phase diagram vanishes. 2A), implying that I-III (line b) must pass through this point; this is not possible with the data from Table 1.

Concluding remarks

Topological phase diagrams do not necessarily provide an accurate description of the position of triple points and of two-phase equilibria. In fact, those positions are as accurate as experimental data and extrapolations can cause them to be. Nonetheless, the phase behavior and in particular the stability hierarchy of the phases involved can be accurately reproduced in a topological phase diagram. This may seem contradictive, but it implies that the overall position of the stability domains of the different phases relative to each other can be accurately determined. To ensure the accurate positioning of stable phase domains, the consistency within a topological phase diagram needs to be verified.

In the present paper, a number of feedback systems for the construction of topological phase diagrams has been highlighted. Thermodynamic cycles for state functions such as

Figure 1 .

 1 Figure 1. Annotated version of the topological pressure-temperature phase diagram presented in Figure 4 of reference [1]. The stable phase regions for forms I, II, and III, are indicated in

  Figure1in reference[1] (see the annotated Figure A.1 in the Appendix of the present paper), it can be seen that the peak is actually of the same magnitude as the solid-solid transitions observed in curve 'a' (III→I) and curve 'e' (II→I) (see Figure A.1 in the Appendix below).

P

  I-liq (/MPa) = 5.57 T I-liq (/K) -2355 (2) P II-liq (/MPa) = 4.38 T II-liq (/K) -1815 (3) P III-liq (/MPa) = 4.29 T III-liq (/K) -1770

Figure 2 .

 2 Figure 2. (A). The topological P-T phase diagram redrawn after the one presented in Figure 3 of reference [1]. The triple points have been numbered as indicated in the original figure. The equilibrium curves can be found through the triple points they cross, 5-9-6: I-II, 4-9-7: I-III, 10-9-8: II-III. The stability hierarchy of the equilibria has not been indicated. (B). The pressure-temperature phase diagram of piracetam with the two-phase equilibria obtained with the data in Table1. The solid-solid equilibria are marked by letters a: I-II (slope 2.08 MPa K - 1 ), b: I-III (1.67 MPa K -1 ), and c: II-III (2.71 MPa K -1 ). In addition, the solid-liquid equilibria I-liq (5.57 MPa K -1 ), II-liq (4.38 MPa K -1 ), and III-liq (4.29 MPa K -1 ) are represented. Point k is the intersection of I-liq and III-liq, which ought to coincide with the triple point I-III-liq (point 7 in Figure2A), implying that I-III (line b) must pass through this point; this is not

*

  enthalpy and specific volume allow for example the determination of missing enthalpy changes between phases as demonstrated above. Moreover, hierarchies within the phase diagram must remain consistent, for example the hierarchy of the slopes of the different transitions dP/dT I-III < dP/dT I-II < dP/dT II-III , which is implied by the monotonicity of the equilibrium curves[2], but as well for experimental data such as T II-III < T I-II < T I-III and ΔH III→II < ΔH II→I < ΔH III→I . The consistency of the measurement data may be obvious, but it does function as a feedback system, when extrapolations lead to other results. What to do when feedback fails to confirm consistency depends on each individual case, but a first reflex should be a thorough review of the experimental data. Since this transition is not observed, its enthalpy can be calculated from experimental data with ∆H III→I + ∆H I→II + ∆H II→III = 0.

Figure A. 1 .

 1 Figure A.1. Annotated differential scanning calorimetry curves for piracetam form III on heating following Figure 1 from reference [1].

Table 1 . Data from the text and from Tables 4 and 7 in reference [1] for the phase transitions of the three most prevalent forms of piracetam Transition T a /°C T b /K ∆H a /J g -1 ∆V f /cm 3 mol -1 ∆V g /cm 3 g -1 (dP/dT) h /MPa K -1

 1

	II→I	91.9 c 87 d	365 360	24.8 -	-4.71	-0.0331	-2.08
	III→I	100.9 c 97 d	374 370	26.8 -	-6.18	-0.0435	-1.67
	III→liq	139.4 c	412.4	211.3	16.97	0.1194	4.29
	II→liq	141.5 e	414.5	211.8 e	16.59	0.1167	4.38
	I→liq	149.8 c	422.8	187	11.29	0.07942	5.57
	a From Table

Table 2 . Coordinates of triple points I-II-liq, I-III-liq, II-III-liq, and I-II-III derived from the data provided in Table 1

 2

	Triple point	P /MPa	T/K	Position in Figures
				2A and 2B
	I-II-liq	170	453	Point 6
	I-III-liq	195	458	Point 7
	II-III-liq	457	518	Point 8
	I-II-III	-84	320	Point 9

Acknowledgements

This work was partially supported by the Spanish Ministry of Science and Innovation through the project FIS2014-54734-P and by the Generalitat de Catalunya under the project 2014 SGR-581.

Appendix

Estimation of the enthalpy changes from the DSC curves reported in Figure 1 of reference [1] by weighing the peaks or measuring their areas as triangles, both results related to the area of peak A (the melting peak of Form I) taken equal to 100% of 187 J/g. The method leads to results similar to those provided in reference [1], even if slightly less accurate.


Dayang Wang

Qian Yin
email: yinqian@bnu.edu.cn

Hongfeng Wang

Weighting Features Before Applying Machine Learning Methods to Pulsar Search

Keywords: Weighting, Machine Learning, Pulsar Search, WEKA

In recent years, different Artificial Intelligence methods have been applied to pulsar search, such as Artificial Neural Network method, PEACE Sorting Algorithm, Real-time Classification method. In this paper, Weighting Feature method before applying machine learning (ML) was proposed. We give weight to each feature according to its ability to distinguish pulsar and non-pulsar candidates. The ability is determined by the separation degree of the distribution of pulsars and non-pulsars on particular feature. And then use the ML methods to classify different types of candidates. The results show that this method is significant. The accuracy of identifying pulsars and modeling time were both improved after weighting.

Introduction

Pulsar is fast rotated neutron star, which periodically sends pulse signal whose period is short and very stable. Pulsar plays an important part in physics, astronomy and many other fields. In recent years, AI methods like image pattern recognition [START_REF] Zhu1 | [END_REF], artificial neural network method and scheduling algorithm are used in pulsar search. Lee et al. (2013) proposed the PEACE sorting algorithm to search pulsar, which had obtained good results. Lyon et al. (2016) used the GH-VFDT (Gauss-Hellinger Very Fast Decision Tree) to distinguish the candidate, with recognition rate of pulsars over 90% [START_REF] Xu | Application of Artificial Intelligence in Pulsar Screening[END_REF]. While GH-VFDT obtained a high recognition rate of pulsars, the difference between the abilities of different features to distinguish the pulse and non-pulsar are not reflected. Thus, in this paper, we add different weights to the eight features before the machine learning process according to their separation degree. Results show that weighting improves both the accuracy rate of classification and modeling time.

The structure of this paper is shown as follows: the related work is mentioned in section 2; the Feature Weighting method is proposed in section 3; and with its corresponding experiments are showed and analyzed in section 4 and 5; the section of conclusion comes as the end.

2

Related Work

Feature

In the process of searching for pulsar signals with radio telescope, the most basic data are obtained. These data are subjected to Removing signal interference, de dispersion, FFT [START_REF] Zhou | X-ray pulsar cycle fast search algorithm for nonequal interval time data[END_REF] [START_REF] Lyon 1⋆ | Fifty Years of Pulsar Candidate Selection: From simple filters to a new principled real-time classification approach[END_REF]. Pulsar Feature Lab and Presto [START_REF] Lee | [END_REF] are used to process the primitive data into these eight features.

Dataset

Three separate datasets were used to the measure the performance of ML methods on pulsar search. The small scale dataset is LOTAAS which was obtained during the LOTAAS survey and is currently private. The medium scale dataset HTRU2 was obtained during an analysis of HTRU Medium Latitude data by Thornton (2013). The large scale dataset HTRU1 is produced by Morello et al. The detailed information of the three datasets is summarized in the table1. It is obvious that when we are classifying a pulsar candidate via its feature, the feature that has a high degree of separation between pulsars and non-pulsars weighs more than other features. Therefore, this paper naturally adds different weights to the eight features according to their separation degree between different types of candidates. As a specific feature, this paper defines the separation degree as follows:

2 P I l l Ab R R     (1)
In this formula, as can be seen from figure2, for a particular feature, Ab denotes the separation degree, l means the coincident area of pulsar and non-pulsar, Rp denotes the width of the distribution of the pulsars on the feature, while R I means the distribution width of non-pulsars. The distribution of features between candidates can be considered as natural distribution. According 3σ principle, features of almost all candidates will be within the range of feature box. By analyzing the data from LOTAAS, HTRU2 and HTRU1, this paper get the weight of each feature W i ( i = 1~8).

W= (1.64, 0.85, 1.91, 1.38, 1.29, 1.50, 2.03, 1.18)

(2)

Experiments

In this part, weighting each feature before utilizing ML methods on the datasets are proposed. Classification accuracy and modelling time are both taken to be criterion to judge the performance of the methods. The paper supposes weighting is useful if methods improves the accuracy or improves the modeling time. What's more, accuracy goes before modeling time. In conclusion, for the five ML methods SMO, IBK, JRIP, J48 and RandomForest, weighting either improves the accuracy or modeling time, or in the worst cases, weighting will at least be the same as not weighting.

Discussions

This part explains why SMO, IBK, JRIP, J48 and RandomForest are selected to test the effects of weighting instead of other ML methods. In this paper, we actually experimented various ML methods using WEKA.

Conclusion

Due to its stable cycle, Pulsar plays a very important part in physics, astronomy and many other fields. Traditional ways of pulsar search are manual. In recent years, Artificial intelligence is widely used in various fields and achieves great success. Therefore, AI methods are gradually applied to pulsar search. This paper is based on the principled real-time classification approach. Eight features are used to describe a pulsar candidate.

Before applying ML methods on datasets, this paper weights each feature according to their separation degree, and then find out that either the accuracy or modeling time is improved after weighting.

Fig. 1 .

 1 Fig. 1. Feature distribution of pulsars and non-pulsars

Fig. 2 .

 2 Fig. 2. Distribution of pulsar and non-pulsars and their coincident area

  and periodic search. Then a pulsar candidate is generated which has some basic Features. Lyon et al. (2016) used eight new features to describe the pulsar candidate. The eight features are Mean of the integrated profile Prof μ , Standard deviation of the integrated profile Prof σ , Excess kurtosis of the integrated profile Prof k , Skewness of the integrated profile Prof s , Mean of the DM-SNR curve DM μ , Standard deviation of the DM-SNR curve DM σ , Excess kurtosis of the DM-SNR curve DM k , Skewness of the DM-SNR curve DM s

Table 1 .

 1 Three pulsar candidate datasetsAnalyzing the statistic distribution of the eight features from the sample data of pulsars and non-pulsars, feature data was extracted from 90, 000 labelled pulsar candidates produced by Morello et al. (2014), via Pulsar Feature Lab. As it is showed in figure1, the data were scaled to the interval of [0, 1]. For each feature, there are two box plots. The orange red box shows the feature distribution of known pulsars, while the box in light blue describes the RFI/noise distribution.

		Dataset	Creator	Time	Volume	Pulsar	Non-pulsar
		LOTAAS Morello et al	2012	5053	66	4987
		HTRU2	Thornton	2013	17898	1639	16259
		HTRU1	Lofar	2013	91191	1196	89995
	3	Methodology-Feature Weighting

Table 2 .

 2 Accuracy and modeling time before and after weighting for LOTAAS.

	Methods	Before	Accuracy	After	Modeling time Before After
	SMO	99.7625%	99.8813%	0.05	0.02
	IBK	99.8417%	99.8615%	0	0
	JRIP	99.8417%	99.8615%	0.08	0.11
	J48	99.8615%	99.8615%	0.05	0.03
	RandomForest 99.8615%	99.8615%	0.7	0.47
	For small scale dataset LOTAAS, the experimental results are shown in Table2, after
	weighting, accuracy rate of SMO, IBK and JRIP are improved. Modeling time of J48
	and RandomForest are improved.

Table 3 .

 3 Accuracy and modeling time before and after weighting for HTRU2. For medium scale dataset HTRU2, the experimental results are shown in Table3, after weighting, accuracy rate of SMO is improved. Modeling time of JRIP, J48 and Random-Forest are improved, while IBK remains the same.

	Methods	Before	Accuracy	After	Modeling time Before After
	SMO	97.5640%	97.5696 %	0.22	0.52
	IBK	97.1449 %	97.1449 %	0	0
	JRIP	97.8154 %	97.8154 %	1.84	1.44
	J48	97.8433 %	97.8433 %	0.28	0.2
	RandomForest	97.9942 %	97.9942 %	9.94	7.16

Table 4 .

 4 Accuracy and modeling time before and after weighting for HTRU1.

	Methods	Before	Accuracy	After	Modeling time Before After
	SMO	99.5866 %	99.5866%	0.36	0.36
	IBK	99.5340 %	99.5340%	0.06	0.01

Table 5 .

 5 Accuracy rate of pulsar recognition of various ML methods before weighting In table5, the purple number means the corresponding methods performs better than others. As is shown, with the scale of datasets becomes larger, SMO, IBK, JRIP, J48 and RandomForest have better performance over other algorithms.

	Types	Methods	LOTAAS	HTRU2	HTRU1
	bayes	NaiveBayes	99.5448%	94.4966 %	98.9155 %
	func-	LibSVM	98.7928%	91.1443 %	98.8585 %
	tions	SMO	99.7625%	97.5640%	99.5866 %
	lazy	IBK LWL	99.8417% 99.9010%	97.1449 % 97.7539 %	99.5340 % 99.5175 %
	meta	AdaBoostML Stacking	99.8615% 98.6938%	97.6534 % 90.8426 %	99.5175 % 98.6885 %
	misc	InputMappedClassifier 98.6938%	90.8426 %	98.6885 %
	rules	JRIP	99.8417%	97.8154 %	99.6129 %
		HoeffdingTree	99.8417%	97.4411 %	99.5822 %
	trees	J48 RandomForest	99.8615% 99.8615%	97.8433 % 97.9942 %	99.6074 % 99.6721 %
		RandomTree	99.8021%	96.8432 %	99.5241 %

Acknowledgments. The research work in this paper was supported by the grants from National Natural Science Foundation of China (61472043, 61375045) and the Joint Research Fund in Astronomy (U1531242) under cooperative agreement between the NSFC and CAS, Beijing Natural Science Foundation (4142030). Prof. Qian Yin is the author to whom all the correspondence should be addressed.


Hang Yu
email: yuhang@mail.bnu.edu.cn

Qian Yin
email: yinqian@bnu.edu.cn

Ping Guo
email: pguo@ieee.org

Enlightening the Relationship between Distribution and Regression Fitting

Keywords: Statistical Distribution, Regression, Model Selection, Goodness of Fit

Statistical distribution fitting and regression fitting are both classic methods to model data. There are slight connections and differences between them, as a result they outperform each other in different cases. A analysis model for processing natural data, say astronomical pulsar data in this paper, is proposed to improve data fitting method performance. Then the insight behind the comprehensive fitting model is given and discussed.

Introduction

Statistics is a subject aiming to model the world by rather collecting and analyzing data than to inference and prove precisely. Lots of natural data have been fitted into statistical distributions with significant performances. And with these statistical distributions scientists can utilize their off-the-shelf distribution properties. While unfortunately the diversity of our world is far more complex than the level those statistical distributions can totally model. Alternatively, polynomial regression is also an approach to automatically find empirical laws from data [START_REF] Edwards | On the use of polynomial regression equations as an alternative to difference scores in organizational research[END_REF]. It is proposed to fit various patterns of data, and has already shown its powerful fitting ability. Researchers who use polynomial regression are conveniently not required to preliminarily estimate data distribution [START_REF] Theil | A rank-invariant method of linear and polynomial regression analysis, 3; confidence regions for the parameters of polynomial regression equations[END_REF]. Regrettably, by implementing polynomial regression we eventually get a polynomial expression with less prior properties than traditional statistical distributions have.

Given that features of natural scientific data are relatively fixed, and they will still be studied in the future, it is significant to determine the best fitting method for certain scientific data feature. And we creatively name it fitting label. As a solution to the dilemma, we propose a discriminative fitting model to help scientists automatically fitting scientific data into the optimal expression. Once scientists get the fitting label of a certain feature of natural scientific data, those later scientists are able to directly update the fitting expression parameters on new data's arrival. We operated experiments on an astronomical data set, ANTF Pulsar Catalogue [START_REF] Manchester R N | The Australia telescope national facility pulsar catalogue[END_REF]. Besides, given the reliability of P value is still controversial [START_REF] Nuzzo | Statistical errors[END_REF], and mean squared error and KL divergence are trustworthy measures for both fields of statistics and machine learning, we complementally took them in consideration.

The structure of this paper is shown as follows: the comprehensive fitting model is proposed in section 2; and with it its corresponding experiments are showed and discussed in section 3 and 4; the section of conclusion comes as the end.

Methodology

Comprehensive Fitting Model

As is shown in Algorithm 1, the selection model consists of mainly 3 steps. Step 1 is to determine whether the data is subject to a certain natural distribution or not through significance testing. Step 2, fitting data into significant distributions and operating regression, then calculating their respective similarities. Step 3 is to classify the data's fitting label with distance-based classifier according to the similarities. Consequently, scientists are able to learn as more properties as possible of the data. 1:

Algorithm 1 Comprehensive Model

Π ← Algorithm2 (D, Π, Γ , Ω) 2: Σ ← Algorithm3 D, Σ, Π 3: Fitting Labels π * ← kN N Σ return π *

Statistical Significance Testing

As is shown in Algorithm 2, firstly we need human experts, astronomers for example, to give us candidate statistical distribution families for each feature. These distribution families with empirical parameters are then be tested one by one for statistical significance under selected statistical significance measures, usually P value. Concretely with specified natural science feature data, we perform statistical significance test to screen out distribution families that have statistically large significance levels for the feature. Candidate statistical families and selected significance measures are given into the testing process, If the output P value reaches the threshold, a significance level of 0.05 by convention, the distribution will be recorded as a significant distribution.

Ψ i ← Ψ i ∪ ψi j ; 5:
end for 6:

ψ i ← Ω × Ψ i ; 7: if ψ i τ then 8: Π ← Π ∪ πi; 9:
end if 10: end for return Π

Similarity Evaluation

For this step, we use similarity measures to detect how well the expression fits with the data, which is shown in Algorithm 3. Significant statistical distributions obtained from last step are introduced here, and with it selected similarity measures. We perform t times k-fold cross-validation test to acquire the average similarity between the distribution and scientific data.

Distance-based Classification

It is the final outcome, the fitting label for the feature data, producing step. We can take those similarity measures as distances, and run a distance-based classifier to eventually select the expression with minimal distance as the optimal fitting label. Here we bring a classic classifier k-Nearest Neighbours(abbreviated as kNN) with its parameter k chosen 1.

Experiments

A logical starting point is to test a distribution's normality for its popularity in natural world [START_REF] Faucher-Giguere C A | Birth and evolution of isolated radio pulsars[END_REF]. While for features whose scatter diagram is shown in Fig. 1 seem obviously subject to lognormal distribution [6][7]. We took the logarithm of that data first, then tested the resulting data's normality, which is shown in Fig. 2. The red curve represented statistical distribution fitting, while the green one for polynomial regression fitting. As we know, logarithm of scale parameter and shape parameter of lognormal distribution are respectively expectation and standard deviation of normal distribution. Significance results are shown in Table 1. For those classic statistical testing methods, the values recorded are their P values. Results shows that lognormal has the highest significance among expert-selecting distributions. Similarity evaluation results are shown in Table 2. The polynomial order is empirically chosen 2. Thus lognormal distribution has the lowest relative entropy and mean squared error. So it is the optimal fitting label for pulsar period feature. The temporary parameters of its corresponding normal distribution are mean -0.53455956344351729 and standard deviation 0.99686283436845236.

Discussions

As we know, the insight is a thought on inductive bias of Occam's razor and the principle of maximum entropy. Occam's razor has been the default heuristic technique in natural science for hundreds of years. It can be interpreted as stating that the one among competing hypotheses with the fewest assumptions should be selected [8][9]. The preference for simplicity also comes to machine learning's avoiding overfitting, a low-order fitting curve is almost always better than a high-order one. On the contrary, the second law of thermodynamics states that the total entropy of an isolated system can only increase over time.

It is seen as an axiom of statistical thermodynamics. Here comes the principle of maximum entropy which states that the probability distribution that best represents the current state of knowledge is the one with largest entropy. and it also emphasizes a natural correspondence between statistical mechanics and information theory [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF]. Numbers of statistical distributions widely used by natural scientists to model scientific data can be derived under the principle of maximum entropy. For any certain mean µ and variance σ 2 , the maximum entropy distribution p x µ, σ 2 is normal distribution N µ, σ 2 [START_REF] Jaynes | Probability theory: The logic of science[END_REF].

Conclusions

For scientists we propose a discriminative fitting model through which diverse scientific data will be fitted into the optimal expression. And with our model as much as possible sample information can be used in future study. For data that is subject to a common statistical distribution, they then will be fitted into that distribution with parameters tuned by numerical optimization methods; comparatively for those who are not we model them by polynomial regression. Among main fields of natural science are still features whose fitting label remains unfixed, so it will be a great work to build a fitting label database with our model for scientists.

Input:

  Scientific Data D; Candidate Statistical Distribution Families Π = {π1, ..., πi, ..., πm}; Selected Significance Testing Methods Γ = {Γ1, ..., Γj, ...Γn}; Selected Similarity Measures Σ = {σ1, ..., σj, ...σn}; Polynomial Regression Order e; Cross-Validation Times t and Fold k. Output: Fitting Labels π * .

Fig. 1 .

 1 Fig. 1. Original Plot

Fig. 2 .

 2 Fig. 2. Logarithm Plot

  Candidate Statistical Distribution Families Π = {π1, ..., πm}; Selected Significance Testing Methods Γ = {Γ1, ..., Γn}; Testing Methods Weights Ω = {ω1, ..., ωn}; Significance Level τ . Output: Significant Statistical Distribution Families Π = {π1, ..., πl }.

	Algorithm 2 Statistical Significance Testing
	1: for each πi ∈ Π do
	2:	for each Γj ∈ Γ do
	3:	ψi

Input: Scientific Data D; j ← Γj (D, πi); 4:

Table 1 .

 1 Significance

	Testing of P0

Table 2 .

 2 Similarity

			Evaluation of P0
	P0	Lognormal	Polynomial
	KL divergence 0.0188633190407 0.565469835801
	MSE	2.49154711009 115.375894222

Acknowledgments. The research work in this paper was supported by the grants from National Natural Science Foundation of China (61472043, 61375045) and the Joint Research Fund in Astronomy (U1531242) under cooperative agreement between the NSFC and CAS, Beijing Natural Science Foundation (4142030). Prof. Qian Yin is the author to whom all the correspondence should be addressed.


Peter G Cook
email: peter.cook@flinders.edu.au

Valentí Rodellas

Aladin Andrisoa

Thomas C Stieglitz

Exchange across the sediment-water interface quantified from porewater radon profiles

Keywords: porewater exchange, recirculated seawater, benthic flux, lake, tracers, radon

Water recirculation through permeable sediments induced by wave action, tidal pumping and currents enhances the exchange of solutes and fine particles between sediments and overlying waters, and can be an important hydro-biogeochemical process. In shallow water, most of the recirculation is likely to be driven by the interaction of wave-driven oscillatory flows with bottom topography which can induce pressure fluctuations at the sediment -water interface on very short timescales.

Tracer-based methods provide the most reliable means for characterizing this short-timescale exchange. However, the commonly applied approaches only provide a direct measure of the tracer flux. Estimating water fluxes requires characterizing the tracer concentration in discharging porewater; this implies collecting porewater samples at shallow depths (usually a few mm, depending on the hydrodynamic dispersivity), which is very difficult with commonly used techniques. In this study, we simulate observed vertical profiles of radon concentration beneath shallow coastal lagoons using a simple water recirculation model that allows us to estimate water exchange fluxes as a function of depth below the sediment-water interface. Estimated water fluxes at the sediment water interface at our site were 0.18 -0.25 m/day, with fluxes decreasing exponentially with depth. Uncertainty in dispersivity is the greatest source of error in exchange flux, and results in an uncertainty of approximately a factor-of-five.

HIGHLIGHTS:

• The radon profile method for estimating porewater exchange is reviewed

• A simple recirculation model is presented to aid in profile interpretation

• Uncertainties in the approach are discussed 1. INTRODUCTION Water recirculation through permeable sediments enhances the exchange of solutes and fine particles between sediments and overlying waters. In particular, it allows for a continuous supply of oxidants and fine particulate and dissolved matter (e.g. dissolved nutrients, bacteria, viruses, phytoplankton) into sediment porewaters, while enhancing the release of degradation products and organisms into overlying waters [START_REF] Huettel | Transport and degradation of phytoplankton in permeable sediment[END_REF][START_REF] Huettel | Flow-induced uptake of particulate matter in permeable sediments[END_REF]. As a consequence, porewater exchange is considered a major contributor to the biogeochemical cycling of surface sediments and overlying waters, particularly in the coastal zone [START_REF] Anschutz | Tidal sands as biogeochemical reactors[END_REF][START_REF] Huettel | Benthic exchange and biogeochemical cycling in permeable sediments[END_REF][START_REF] Huettel | Advective Transport Affecting Metal and Nutrient Distributions and Interfacial Fluxes in Permeable Sediments[END_REF][START_REF] Jahnke | Organic matter remineralization and porewater exchange rates in permeable South Atlantic Bight continental shelf sediments[END_REF]. Porewater exchange increases in importance in highly-permeable sandy sediments, which cover >40% of coastal and shelf areas worldwide [START_REF] Riedl | The subtidal pump: a mechanism of interstitial water exchange by wave action[END_REF], where this advective transport of solutes can exceed fluxes driven by molecular diffusion by several orders of magnitude [START_REF] Huettel | Porewater flow in permeable sediments[END_REF]. This advective exchange between porewaters and overlying waters is caused by pressure gradients at the sediment-water interface, which might be forced by several mechanisms spanning a range of spatial and temporal scales, including wave and tidal pumping, interaction of bottom currents and seafloor topography, density instabilities or pumping activities of benthic fauna (bio-irrigation, [START_REF] Huettel | Benthic exchange and biogeochemical cycling in permeable sediments[END_REF][START_REF] Santos | The driving forces of porewater and groundwater flow in permeable coastal sediments: A review[END_REF]. For instance, the passage of waves can produce oscillatory flows that interact with bottom topography (e.g. ripples), producing local increases of pressure that drive fluid exchange across the sediment-water interface. In addition, the passage of wave crests and troughs creates pressure gradients over the seafloor that also enhance porewater exchange [START_REF] Riedl | The subtidal pump: a mechanism of interstitial water exchange by wave action[END_REF]Rutgers van der Loeff, 1981;[START_REF] Webster | Wave Enhancement of Diffusivities within Surficial Sediments[END_REF].

Despite the importance of porewater exchange in coastal biogeochemical cycles, it is still not easy to quantify the advective flux of water and solutes in permeable sediments [START_REF] Boudreau | Permeable marine sediments: Overturning an old paradigm[END_REF][START_REF] Rocha | Sandy sediments as active biogeochemical reactors: compound cycling in the fast lane[END_REF]. Common methods to estimate the rate of porewater exchange across the sediment-water interface in permeable sediments include (1) deploying automated seepage meters to monitor the porewater flow into overlying waters (e.g. [START_REF] Jahnke | Benthic flux of biogenic elements on the Southeastern US continental shelf: influence of pore water advective transport and benthic microalgae[END_REF][START_REF] Cable | Exonerating Bernoulli? On evaluating the physical and biological processes affecting marine seepage meter measurements[END_REF],

(2) constructing mass balances in overlying waters for a tracer supplied by porewater inputs (e.g. [START_REF] Stieglitz | Karstic groundwater discharge and seawater recirculation through sediments in shallow coastal Mediterranean lagoons, determined from water, salt and radon budgets[END_REF]Rodellas et al., 2017), (3) injecting artificial tracers or dye into or above the sediments to trace fluid advection across the sediment-water interface (e.g. [START_REF] Reimers | In situ measurements of advective solute transport in permeable shelf sands[END_REF]Precht and Huettel, 2004), and (4) modeling the depth profiles of temperature, electrical conductivity or dissolved species in sediments to evaluate the sediment-water exchange rates (e.g. [START_REF] Cable | In situ evaluation of nearshore marine and fresh pore water transport into Flamengo Bay, Brazil[END_REF][START_REF] Savidge | Using a Thermal Proxy to Examine Sediment-Water Exchange in Mid-Continental Shelf Sandy Sediments[END_REF]. Among these dissolved compounds, short-lived, naturally-occurring radionuclides (e.g. 224 Ra and 222 Rn) have been widely used, mainly because their half-lives are well suited to the common timescales of these porewater exchange processes (e.g. [START_REF] Colbert | Shoreline and seafloor fluxes of water and short-lived Ra isotopes to surface water of San Pedro Bay, CA[END_REF][START_REF] Cable | In situ evaluation of nearshore marine and fresh pore water transport into Flamengo Bay, Brazil[END_REF][START_REF] Cai | 224Ra:228Th disequilibrium in coastal sediments: Implications for solute transfer across the sediment-water interface[END_REF]. Radium-224 (half-life = 3.6 d) is produced by radioactive decay of 228 Th, but partitions to aquifer solids which reduces the sensitivity of dissolved 224 Ra activities to porewater exchange fluxes.

As a noble gas, radon ( 222 Rn; half-life = 3.8 d) is an excellent tracer of porewater exchange because it is not affected by chemical and biological processes occurring within sediments. In this study, we review the application of radon porewater profiles in sediments to estimate the porewater exchange across the sedimentwater interface driven by pressure fluctuations on short time-scales, such as those produced by wave action. We simulate observed radon profiles using a numerical model of water recirculation that allows us to estimate porewater exchange rates as a function of depth below the sediment-water interface, and to explore sensitivity of estimated fluxes to model parameters. An alternative simpler approach is also presented, where water fluxes are estimated based on changes in the observed radon concentration gradient with depth.

THEORY

Radon-222 is a natural environmental tracer that has been used for quantifying groundwater inflows to streams and estuaries (e.g. [START_REF] Cook | Quantifying groundwater discharge to Cockburn River, southeastern Australia, using dissolved gas tracers 222 Rn and SF 6[END_REF][START_REF] Genereux | Use of radon-222 and calcium as tracers in a three-end-member mixing model for streamflow generation on the West Fork of Walker Branch Watershed[END_REF]) and the ocean (e.g. [START_REF] Cable | Application of 222Rn and CH4 for assessment of groundwater discharge to the coastal ocean[END_REF]) for almost three decades. Radon is produced in the sediments by the radioactive decay of 226 Ra, which is part of the 238 U decay chain and it is found in the sediment solids and in porewater. After porewater containing radon discharges to surface water bodies, radon activities in the surface water decrease due to gaseous exchange with the atmosphere and radioactive decay. Radon concentrations in surface water are therefore always much lower than concentrations within porewaters. When surface water infiltrates, its dissolved radon concentration will increase, according to

!! !! = γ -λc (1)
where c is the radon concentration, γ is the production rate of dissolved radon, λ is the radon radioactive decay coefficient (0.1818 d -1 ), and t is time. After a residence time of a few weeks the concentration will reach secular equilibrium, in which the rate of production is exactly balanced by the rate of radioactive decay. The concentration at secular equilibrium is equal to γ/λ. If the radon production rate within the sediments is constant, then the radon concentration in porewater beneath the seafloor will increase with depth, up to a depth where secular equilibrium is reached. The depletion of radon at shallow depths may therefore be used to derive the exchange between porewater and overlying waters.

Radon deficit model

The most commonly applied approach to estimate the flux of radon across the sediment-water interface is based on the deficit of porewater radon relative to radon concentrations at secular equilibrium, i.e. the radon concentration that would occur without solute exchange (e.g. [START_REF] Martin | Magnitudes of submarine groundwater discharge from marine and terrestrial sources: Indian River Lagoon, Florida[END_REF][START_REF] Cable | In situ evaluation of nearshore marine and fresh pore water transport into Flamengo Bay, Brazil[END_REF]. This deficit must be equal to the total net flux of radon into overlying waters. More generally, and assuming 1D flow, the net radon flux at depth z' can be written as:

𝐽 !" = 𝛾 -𝜆𝑐 𝜃𝑑𝑧 ! !! (2)
where θ is sediment porosity and c is the radon concentration at each depth z.

Integrating from z'=0 to ∞ gives the net radon flux across the sediment-water interface.

Aside from the sediment porosity and the radon concentrations at the different depths, this approach requires knowledge of the production rate in sediments (γ), which can either be calculated from slurry-equilibration experiments [START_REF] Colbert | Shoreline and seafloor fluxes of water and short-lived Ra isotopes to surface water of San Pedro Bay, CA[END_REF] or derived from deep porewater radon concentrations, which are assumed to be unaffected by porewater exchange [START_REF] Cable | In situ evaluation of nearshore marine and fresh pore water transport into Flamengo Bay, Brazil[END_REF].

Advection cycling model

Exchange of solutes between rivers and lakes and underlying porewaters has frequently been simulated using mass balance models that represent the porewater zone as a perfectly mixed reservoir of constant depth [START_REF] Bencala | Simulation of solute transport in a mountain pool-and-riffle stream with a kinetic mass transfer model for sorption[END_REF][START_REF] Cook | Quantifying groundwater discharge to Cockburn River, southeastern Australia, using dissolved gas tracers 222 Rn and SF 6[END_REF][START_REF] Gooseff | Comparing transient storage modeling and residence time distribution (RTD) analysis in geomorphically varied reaches in the Lookout Creek basin, Oregon, USA[END_REF][START_REF] Stieglitz | Karstic groundwater discharge and seawater recirculation through sediments in shallow coastal Mediterranean lagoons, determined from water, salt and radon budgets[END_REF]. Multiple reservoirs sometimes have been used, although such reservoirs usually operate in parallel so that each reservoir is directly connected to the surface water [START_REF] Choi | Characterizing multiple timescales of stream and storage zone interaction that affect solute fate and transport in streams[END_REF]. In most cases, the focus of these studies has been on reproducing exchange fluxes rather than porewater concentrations, although [START_REF] Lamontagne | Estimation of hyporheic water residence time in situ using 222 Rn disequilibrium[END_REF] used a single reservoir model to relate water fluxes to radon concentrations within porewater of a river hyporheic zone that was assumed to be perfectly mixed.

As the objective of the present study is to interpret vertical radon profiles rather than mean porewater concentrations, we consider a series of vertically-stacked reservoirs rather than a single reservoir. This type of model, sometimes referred to as a compartmental mixing model, has been widely applied in soil water and groundwater transport simulations [START_REF] References Adar | Estimation of spatial recharge distribution using environmental isotopes and hydrochemical data, I. Mathematical model and application to synthetic data[END_REF][START_REF] Harrington | A compartmental mixing-cell approach for the quantitative assessment of groundwater dynamics in the Otway Basin, South Australia[END_REF][START_REF] Kirk | A deuterium-calibrated groundwater flow model of a regional carbonate-alluvial system[END_REF], but not previously in the context of porewater exchange.

Water and solutes are assumed to be perfectly mixed within each cell, and only advective fluxes between the cells are considered. Thus hydrodynamic dispersion is simulated implicitly through mixing within the cells rather than explicitly included in the governing equations. Compartmental mixing cell models have been shown to produce similar results to advective-dispersive models provided that advection is the dominant transport process and that the size of the mixing cell is appropriately chosen [START_REF] Xu | Criteria for resolution-scales and parameterisation of compartmental models of hydrological and ecological mass flows[END_REF].

The model presented here represents a two-dimensional recirculation cell, in which flow reverses periodically. This type of flow system might be produced by the passage of waves across the water surface, with surface water moving into the sediments beneath the wave peaks (high pressure), and exiting beneath the wave troughs (low pressure). We assume that the vertical scale of the recirculation cell is much greater than the horizontal scale, so that horizontal travel times are negligible (this is discussed further below). The two-dimensional model thus collapses into two one-dimensional profiles that exchange water and solutes. The fully saturated porewater zone is discretized into a number of layers (cells), each of which is assumed to be perfectly mixed. Each cell is assumed to continually exchange water with the cells immediately above and below it (the uppermost cell representing the surface water), and also to exchange water with the cell at the equivalent depth in the adjacent profile (Figure 1). The downward flow in one profile is thus balanced by the upward flow in the second profile, and the change in vertical water flux with depth determines the exchange flux between the two profiles. For each profile, we solve

! !" !" = - ! !" 𝑞 ! 𝑐 + 𝛾𝜃 -𝜆𝜃𝑐 -𝑆 (3)
where 𝑞 ! is vertical water flux (which is a function of depth), and S is a mass flux term that allows for flow between the upwelling and downwelling profiles. 𝑞 ! is defined so that downward fluxes are positive and upward fluxes are negative.

Dispersion is not explicitly simulated, but it is implicitly simulated based on the size of the mixing cells. The relationship between cell size (Δz) and implicit dispersivity (α) is given by [START_REF] Kirchner | Applicability of compartmental models for simulating the transport of radionuclides in soil[END_REF][START_REF] Shanahan | Transport in Lake Water Quality Modeling[END_REF]:

𝛼 = ∆! ! (4)
The water flux decreases with depth in the downwelling profile, so that for the first phase of the recirculation cycle, water moves downward on the left hand side and upwards on the right hand side, as shown in Figure 1. Horizontal flows occur from left to right, and are given by

𝑞 ! 𝑖 = 𝑞 ! 𝑖 -1 -𝑞 ! 𝑖 (5)
where 𝑞 ! (𝑖) is the horizontal water flux at the depth represented by cell i, qv(i-1) is the downward water flux into cell i from the overlying cell and qv(i) is the vertical water flux from cell i into the underlying cell. Note that both qv and qh are volumetric fluxes per square area of the sediment surface, and so qh is not a Darcy velocity in the traditional sense. Note also that the term S in Equation 3 can be obtained by multiplying qh(i) by the concentration in cell i of the downwelling profile. S will therefore be positive for downwelling profiles (-S is negative, representing a mass loss) and negative for upwelling profiles (-S is positive, representing a mass gain).

The solute mass balance equations are:

∆𝑐 ! 𝑖 = ∆𝑡 ! ! ! !!! ! ! ! !!! !! ! ! ∆! + 𝛾 -𝜆𝑐 ! 𝑖 (6) ∆𝑐 ! 𝑖 = ∆𝑡 !! ! ! ! ! ! ! !!! !! ! ! ∆! + ! ! ! ! !! ! ! !!! ! ! ! ! !! ! ! ∆! + 𝛾 -𝜆𝑐 ! 𝑖 ( 7
)
where The upper boundary condition is constant concentration (c=c0). The lower boundary is the concentration of radon in equilibrium with the rate of production (c=γ/λ), although in practice the lower boundary is set to be sufficiently deep so that it does not affect simulation results. The key parameters in the model are the concentration in the overlying surface water (c0), the sediment characteristics (θ, γ), and the characteristics of the recirculation, which include the time for a completed cycle (tr), and the velocity profile qv(i).

It should be noted that as the period of the recirculation (tr) becomes small, radon concentrations derived from the oscillating flow model are equivalent to those observed in a simple 1D steady state model with flow occurring simultaneously in both directions. This is given by

! ! !!! ! ! !!! !! ! ∆! - ! ! ! ! ! ! !! !!! ∆! + 𝛾 -𝜆𝑐 = 0 (8)

Dispersion model

For small recirculation times (tr), the solute profiles in a reversing flow field can also be expressed in terms of an enhanced dispersion coefficient, rather than by directly simulating advection [START_REF] Qian | Depth-Dependent Dispersion Coefficient for Modeling of Vertical Solute Exchange in a Lake Bed under Surface Waves[END_REF]. The flux of radon at any depth can be expressed

𝐽 !" = 𝐷 ! 𝜃 !" !" ( 9
)
where JRn is the radon flux and De is the enhanced dispersion coefficient (m 2 d -1 ).

Qian et al. ( 2009) developed a two-dimensional hydraulic model to examine the effect of wave action on porewater solute profiles, and showed that the value of the enhanced dispersion coefficient decreased exponentially with depth, and could be approximated by

𝐷 ! (𝑧) = !!"# !" 𝑒𝑥𝑝 !!.!"! ! ( 10
)
where K is the sediment hydraulic conductivity, a is the half-wave amplitude and L is the wavelength (m).

METHODS

3.1. Radon sampling: study site, sampling and analyses

Porewater profiles for radon analysis were collected at La Palme lagoon, located on the western French Mediterranean coastline. La Palme is a small (5 km 2 surface area), shallow lagoon, with mean and maximum water depths of 0.6 and 1.5 m respectively. It is connected with the Mediterranean Sea through a small opening in the coastal sand spit and it receives continuous fresh groundwater inputs mainly from a regional karst aquifer [START_REF] Stieglitz | Karstic groundwater discharge and seawater recirculation through sediments in shallow coastal Mediterranean lagoons, determined from water, salt and radon budgets[END_REF]. The internal mixing of the lagoon and its exchange with coastal waters is driven primarily by the strong north-westerly winds characteristic of the region (regularly exceeding 10 m/s).

Tidal forcing plays a minor role on the hydrodynamic functioning of this lagoon (tidal variations in the Mediterranean Sea are usually less than 0.4 m and exchange between the lagoon and the sea is highly restricted). Most of the lagoon is covered by coarse-grained highly-permeable sediments. A recent study conducted by [START_REF] Stieglitz | Karstic groundwater discharge and seawater recirculation through sediments in shallow coastal Mediterranean lagoons, determined from water, salt and radon budgets[END_REF] revealed that wind-driven horizontal pressure gradients at the sediment-water interface produce the recirculation of large amounts of lagoon water through surface sediments. Indeed, they estimated that the equivalent of the volume of the entire lagoon recirculates through the sediments every 25 days. La

Palme is thus an ideal site to evaluate the exchange of porewater across the sediment-water interface by using radon porewater profiles.

Porewater samples for radon analysis were collected at 2 different locations (Pz1 at 42.9741°S, 3.0163°E and Pz2 at 42.9391°S, 3.0248°E) using a drive point piezometer. At each location, samples were collected at depths of 0.05 (only Pz1), 0.10, 0.15, 0.20, 0.30, 0.50, 0.80 and 1.30 (only Pz2) m below the sediment -water interface. 10-mL porewater samples were collected using a gas-tight syringe coupled to the piezometer tubing (minimizing water-air contact) and transferred to 20-mL vials prefilled with a 10-mL high-efficiency mineral oil scintillation cocktail [START_REF] Cable | In situ evaluation of nearshore marine and fresh pore water transport into Flamengo Bay, Brazil[END_REF]. Concentrations of radon were analyzed by liquid scintillation counting on a Quantulus 1220 with alpha-beta discrimination counting (background of 0.2-0.4 cpm; efficiency of 1.6-2.2, depending on the quenching factor of the sample). Surface water samples (2 L) were collected using a small submersible pump to minimize the gas loss and analyzed using the radonin-air monitor RAD7 coupled to an extraction system. Samples were decaycorrected to the time of collection.

Model

The advective compartmental mixing cell model was programmed into Fortran 95.

We used a cell size of Δz = 0.01 m, which is equivalent to an implicit dispersivity of α = 0.005 m. The latter is consistent with a flow length of approximately 0.5 m [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF]. (This is the approximate depth of radon depletion apparent in the measured profiles, and hence also the apparent depth of recirculation.)

Temporal discretization was 10 -7 days (0.0086 seconds), but identical results were obtained using smaller values. The model was run for at least 20 days, so that radon concentrations reach dynamic equilibrium and are unaffected by initial conditions. Sediment porosity and radon production rate were both assumed to be constant (i.e., do not vary with depth).

RESULTS

The two porewater profiles Pz1 and Pz2 showed radon concentrations increasing rapidly with depth up to around 30 -50 cm depth. Radon concentrations below these depths were relatively constant, with maximum concentrations of approximately 5000 Bq m -3 and 2500 Bq m -3 for Pz1 and Pz2 respectively, which likely reflect concentrations reaching secular equilibrium. The lower value observed at 0.8 m for PZ1 is likely due to analytical or sampling uncertainty or may reflect a spatial variation in the radon production rate. The difference between the measured radon concentrations at shallow depth and these equilibrium concentrations sustained by the production rate indicates that there is a significant exchange of radon between porewaters and overlying waters.

Deficit model

The radon production rates for the two sites are approximately 900 and 450 Bq m 3 d -1 for Pz1 and Pz2, respectively, as derived from the maximum concentrations measured at each site (assuming constant production rates over depth). By applying equation 2 and using a porosity (θ) of 0.4, we estimated a total net flux of radon into overlying waters of 58 and 38 Bq m -2 d -1 for Pz1 and Pz2, respectively.

The radon flux estimated here refers to the total net loss of radon from sediments into surface waters and thus includes also the flux of radon supplied by molecular diffusion. The net advective-dispersive flux of radon from sediments can be calculated as the difference between the estimated total flux and the diffusive flux, which can be approximated using Fick's First Law and radon diffusion coefficients corrected for both temperature and tortuosity (~1 × 10 -4 m 2 d -1 ; [START_REF] Boudreau | Diagenetic models and their implementation: modelling transport and reactions in aquatic sediments[END_REF][START_REF] Li | Diffusion of ions in sea water and in deep-sea sediments[END_REF]). By using the radon gradients over depth measured in the shallowest porewaters (above 20 cm), the diffusive flux is on the order of 10 -1 Bq m -2 d -1 , which is 2 orders of magnitude lower than the total radon flux. Therefore, radon diffusive fluxes in the studied profiles can be neglected and the total 222 Rn flux can be attributed to advective-dispersive fluxes.

Advection cycling model

Observed radon profiles were simulated using constant surface water concentrations of 30 (Pz1) and 90 (Pz2) Bq m -3 , and 𝑞 ! 𝑧 and 𝑡 ! were varied in a trial-and-error fashion, until good fits with observed profiles were obtained. It was found that best-fits to radon profiles were produced with water fluxes that decreased exponentially with depth, and so this was adopted in all simulations.

The simulated radon profile for Pz1 is shown in Figure 3, and closely matches the observed profile to a depth of 0.5 m. This simulation uses a recirculation time of tr/2 = 10 -5 days (0.864 seconds), although identical profiles are produced for recirculation times on the order of <10 -2 days. For longer times, different upwelling and downwelling profiles are obtained (Figure 4). In this case, the profile that has just completed its downwelling phase (orange line in Figure 4) has much lower concentrations than the upwelling profile (blue line). This essentially represents differences that would be observed depending on the sampling time in relation to the phase of the cycling. values (0.15 -0.5 m) are towards the lower limit. The reverse applies for the higher values of k. This is extremely unlikely to be the case (probability of 2 -6 ; or less than 2%). Thus, the true uncertainty in k is much less than this, and probably closer to 2.0 < k < 10 m -1 . Values of k are most sensitive to radon concentrations at 0.2 -0.4 m depth (Figure 5), and increasing the precision of these measurements would significantly improve the accuracy of the flux attenuation rate estimate.

Radon concentrations in profile Pz2 are reproduced using a production rate of γ = 450 Bq m -3 d -1 , and surface flux 𝑞 ! 0 = 0.18, but with other parameters identical (θ = 0.4; 𝑡 ! /2= 10 -5 d (0.86 s); with k = 5 m -1 ).

Dispersion Model

It is also possible to represent the observed radon data in terms of enhanced dispersion coefficients rather than explicitly considering advective fluxes. By equating Equations 2 and 9, we obtain

𝐷 ! 𝜃 !" !" !!!! = γ -𝜆𝑐 ! !! 𝜃𝑑𝑧 (11)
The dispersivity profile can therefore be obtained by approximating !"

!" using the observed changes in radon concentration between sampling depths and solving Equation 11.

The profile of enhanced dispersion coefficient calculated for profile Pz1 in this way is depicted in Figure 7a. The dispersion coefficient is related to the flux by 𝐷 ! = 𝑣 𝛼, where ν is the water velocity. As dispersion will occur under both upward and downward flow, it is therefore related to the recirculation flux according to

𝐷 ! = ! ! ! ! = !! ! ! ! (12)
where 𝑞 ! is the upward or downward flux that occurs during respective upwelling and downwelling phases, and 𝑞 ! is the mean upward or downward flux averaged across the two phases. Thus we can calculate the average flux 𝑞 ! by rearranging equation 12. This is shown in Figure 7b, where it is also compared with the flux profile obtained from the advection cycling model. The surface flux calculated using Equations 11 and 12 is 0.31 m d -1 , which is similar to that calculated from the mixing cell model (0.25 m d -1 ).

DISCUSSION

Model sensitivity

The results from both the advection cycling and the dispersion model are a function of the selected longitudinal dispersivity (expressed as cell size for the advection cycling model). Changing the selected dispersivity will result in proportional changes on the estimated water fluxes. Since dispersivity is a scale dependent parameter, longitudinal dispersivity (α) has been often related to the length of the flow path at the field scale [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF][START_REF] Neuman | Universal scaling of hydraulic conductivities and dispersivities in geologic media[END_REF].

Assuming a flow path length of 0.5 m, we would obtain an approximate dispersivity of 0.005 m, which is the value used in this study. The estimated dispersivity would largely depend on the length scale selected, which may be difficult to define when dealing with porewater fluxes. As an alternative approximation, [START_REF] Qian | A vertical dispersion model for solute exchange induced by underflow and periodic hyporheic flow in a stream gravel bed[END_REF] suggested that pore-scale dispersivities in a sediment bed can be approximated by the average particle size. Since the sediment of the study site are mainly composed by coarse-grained particles, the mean particle size is likely ranging from half to a few millimeters, and the expected range of α would therefore be on the order of 0.001 m. This gives a factor-of-five uncertainty in α, and hence also in water flux.

Aside from longitudinal dispersivity, which is the single most important control on model results, the advection cycling model is also sensitive to other input parameters. Although exchange flux (𝑞 ! ) is proportional to porosity (θ), this parameter is relatively easy to estimate in unconsolidated sediments, and it is not expected to vary in a wide range. The radon production rate (γ) is constrained by the radon concentrations measured in deep porewaters, and so it is important for samples to be collected from sufficient depth. In our example, γ was assumed to be constant with depth, although variations in θ and γ with depth can be introduced to the model if they are expected to occur. In this case, γ would need to be derived from equilibration experiments, and previous studies have shown that these measurements can have significant uncertainties [START_REF] Key | Emanation of radon-222 from marine sediments[END_REF][START_REF] Berelson | Radon-222 as a tracer for mixing in the water column and benthic exchange in the southern California borderland[END_REF]. Model-derived water fluxes in the uppermost layers of the sediment are relatively insensitive to the selected production rate (γ), but γ exerts an important control on the water fluxes simulated for deeper layers, where radon concentrations approach secular equilibrium. Accurate estimation of water fluxes at depth, and the rate (and functional form) of flux attenuation, is also highly sensitive to the analytical precision of radon measurements. We assumed that 𝑞 ! decreased exponentially with depth, as this is commonly assumed (e.g. [START_REF] Qian | Depth-Dependent Dispersion Coefficient for Modeling of Vertical Solute Exchange in a Lake Bed under Surface Waves[END_REF][START_REF] Fram | Miniature thermistor chain for determining surficial sediment porewater advection[END_REF][START_REF] Wilson | Using heat as a tracer to estimate the depth of rapid porewater advection below the sediment-water interface[END_REF] and because the exponential decrease produced the best fit to the data. Finally, the model is not very sensitive to variations in the time for each phase for the cycle when tr is on the order of <10 -2 days. For a given value of dispersivity, the uncertainty in all these input parameters and analytical precision result in an uncertainty in the water flux across the sediment -water interface (𝑞 ! 0 ) of approximately a factor-of-two.

Thus, as previously discussed, uncertainty in dispersivity dominates the uncertainty in surface water flux.

A limitation of the model is that it neglects horizontal travel times between downwelling and upwelling profiles. This may be reasonable for small-scale recirculation driven by wave action, as at any time upwelling and downwelling zones are likely to be separated by only a few centimetres. The assumption is less likely to be valid for recirculation systems driven by processes operating over longer timescales, such as seiches or tides. For these, porewater concentrations (such as those shown in Figure 4) are likely to be underestimated, particularly for upwelling profiles.

Quantification of water exchange across the sediment-water interface

A number of previous studies have estimated net radon fluxes across the sediment-water interface by integrating the radon deficit in sediment porewater (Equation 2) (e.g. [START_REF] Cable | In situ evaluation of nearshore marine and fresh pore water transport into Flamengo Bay, Brazil[END_REF][START_REF] Martin | Magnitudes of submarine groundwater discharge from marine and terrestrial sources: Indian River Lagoon, Florida[END_REF]. However, it is not straightforward to calculate the water flow from this data. The common approach to convert this radon flux into specific porewater discharge is by dividing the mass flux by the radon concentration in the shallowest porewater sample [START_REF] Cable | In situ evaluation of nearshore marine and fresh pore water transport into Flamengo Bay, Brazil[END_REF][START_REF] Martin | Magnitudes of submarine groundwater discharge from marine and terrestrial sources: Indian River Lagoon, Florida[END_REF]. Given that it is extremely difficult to collect porewater samples for radon analysis in the top centimeters of the sediment, the nearest-surface porewater samples are commonly collected at 0.05-0.1 m depth. In the case of the profiles collected in La Palme lagoon, the radon concentration in the uppermost porewater sample in PZ1, which was collected at a depth of 0.05 m, is 930 Bq m -3 . Had we followed this approach to estimate the water exchange between porewater and overlying waters, we would have obtained a flow of 0.06 m d -1 , which would underestimate the water flow derived from the advection cycling model (𝑞 ! 0 = 0.25 m d -1 ) by a factor of approximately 4 (or by a factor of 8 if the shallowest sample was collected at 0.1 m depth). The main difference between these estimates is that the midpoint of the uppermost cell for the advection cycling model is at a depth of 0.005 m (with a simulated radon concentration of 240 Bq m -3 ), which is a depth virtually impossible to sample for radon analysis. Although the concentration difference between the uppermost cell and the surface water must be considered to estimate the water exchange rate, the surface water concentrations are often small when compared with porewater concentrations. Thus, whilst the commonly applied radon deficit approach allows estimating net radon fluxes across the sediment-water interface, water fluxes are more accurately quantified by modeling the radon distribution with depth. An alternative approach is to calculate the dispersion coefficient by dividing the net radon flux at each depth by the concentration gradient. Even if samples are collected from below 5 cm, the dispersion coefficient calculated in this way does not appear to greatly underestimate the dispersion coefficient at the sedimentwater interface. This is because

!"

!" varies more slowly with depth than does c. The water flux can then be estimated by dividing the surface dispersion coefficient by the dispersivity.

One of the advantages of the radon approach described here is that it allows the estimation of water fluxes as a function of depth. However, as discussed above, an accurate estimation of exchange fluxes at deeper depths would require more precise radon measurements (e.g. counting deep radon porewater samples for longer times), increased sampling resolution and/or radon equilibration experiments [START_REF] Colbert | Temporal and spatial variability of radium in the coastal ocean and its impact on computation of nearshore cross-shelf mixing rates[END_REF] to provide independent estimates of the radon production rate. While most of the porewater exchange studies have focused on the fluxes across the sediment-water interface, determining the exchange fluxes at deeper depths is important for understanding the biogeochemical cycles in sediments. The penetration depth of reactants (e.g. oxygen), for example, will depend on the advective porewater velocities, as well as on the consumption/production rates in the sediment layers (Precht et al., 2004).

The approaches described in this paper are most appropriate in those systems where the driving force generating horizontal pressure gradients at the sedimentwater interface oscillates in relatively short temporal scales (seconds to hours).

Larger recirculation times (hours to days) would result in profiles that would significantly change depending on the sampling time in relation to the phase of the advection cycle (upwelling or downwelling) (Figure 4). The proposed approaches are thus best suited to quantifying porewater exchange fluxes produced by the undulating pressure at the seafloor generated by gravity waves interacting with relatively flat sediment surfaces. These models implicitly include the effects of interaction between wave-driven oscillatory flows and seabed morphology, which may significantly enhance water recirculation through sediments, particularly in areas with a water depth shallower than half the wavelength of the wave [START_REF] Precht | Advective pore-water exchange driven by surface gravity waves and its ecological implications[END_REF]. However, if bedforms (e.g. ripples) are stable on timescales of hours or longer, this might give rise to stable zones of up-and downwelling, and so profiles would vary depending on the area of collection. Note that bottom topography can change significantly over short time scales (e.g. ripple migration), particularly during strong periodic events (e.g. storms) or in areas affected by strong bottom currents (Precht et al., 2004;[START_REF] Savidge | Forcing and Dynamics of Seafloor-Water Column Exchange on a Broad Continental Shelf[END_REF]. Therefore, zones of upwelling and downwelling porewater in permeable sediments would also propagate with ripple migration (Precht et al., 2004). This intensive lateral shifting of up-and downwelling zones within the sediment together with horizontal diffusion and dispersivity may contribute to homogenizing the vertical profiles. In a similar manner, areas of preferential resuspension or deposition of sediments could also release or trap significant volumes of porewaters [START_REF] Santos | The driving forces of porewater and groundwater flow in permeable coastal sediments: A review[END_REF], and thus would also result in significantly different porewater profiles depending on the sampling area. Collecting different radon porewater profiles in the same area should provide additional information on the temporal and spatial scales of the driving forces, by identifying the stability of upwelling and downwelling zones.

The advection cycling and the dispersion models represent thus a reliable method to characterize water exchange across the sediment-water interface driven by pressure gradients reversing at short temporal scales. Radon has advantages over other porewater tracer approaches, as it is more sensitive to low exchange fluxes than temperature [START_REF] Briggs | Practical limits on the use of diurnal temperature signals to quantify groundwater upwelling[END_REF], and is simpler than dye injection approaches. Other methods commonly applied to quantify porewater exchange are not well suited to the estimation of fluxes with such short residence times. In situ seepage meters may alter fluxes above and below the sediment interface due to the presence of the instrument. This might not be significant for fluxes driven by longer-scale pressure changes, but is likely to be important for processes operating on very short time-scales. Tracer mass balances in overlying waters require estimation of the concentration of exchanging water to convert the tracer mass balance to a water mass balance. The appropriate end-member concentration will depend on the hydrodynamic dispersivity for this method, as it does for porewater tracer methods. However, mass balances in overlying waters will have additional uncertainties due to the need to also define other components of the mass balance.

Model-derived results for La Palme Lagoon

The shape of the radon porewater profiles collected in La Palme Lagoon (Figure 2)

suggests that porewater exchange at the sites sampled is driven by pressure gradients reversing at short temporal scales (up to hours). Larger reversing scales would have produced radon profiles similar to those shown in Figure 4. The enhanced diffusion coefficients at the surface water -sediment interface calculated from the advection cycling model (𝐷 ! = 2𝑞 ! 𝛼𝜃 !! = 0.006 m 2 d -1 for Pz1 and 0.004 for Pz2) or estimated from Equation 11(0.008 m 2 d -1 for Pz1 and 0.008 m 2 d -1 for Pz2) can be compared with the relationship derived by [START_REF] Qian | Depth-Dependent Dispersion Coefficient for Modeling of Vertical Solute Exchange in a Lake Bed under Surface Waves[END_REF] and reproduced in Equation 10. Using the simulated attenuation rate (5 m -1 ), we obtain 6.15/L = 5 m -1 , and hence L ≈ 1.2 m. Similarly, !!"# !" ≈ 0.004 -0.008 m 2 d -1 (13)

Using α = 0.005 m, θ = 0.4, and assuming a = 0.05 m then gives K ≈ 1.5 -3.1 m d -1 .

These are typical values of hydraulic conductivity for saturated silty sands [START_REF] Freeze | Groundwater[END_REF], such as those found in La Palme Lagoon. Therefore, the observed surface dispersion coefficient and velocity and their rate of attenuation with depth is consistent with wave action being the principle driver for recirculation in La Palme Lagoon. Since wave dynamics in La Palme Lagoon are mainly controlled by the wind regime, wind forcing appears to be an important driver of porewater exchange in this lagoon, as already suggested by [START_REF] Stieglitz | Karstic groundwater discharge and seawater recirculation through sediments in shallow coastal Mediterranean lagoons, determined from water, salt and radon budgets[END_REF]. Future studies should focus on evaluating the temporal variability of water and solute fluxes across the sediment-water interface in different wind (wave) conditions.

Considering that the water depths of La Palme lagoon usually range from 0.3 to 1.5 m, the figures estimated in this study would imply that the entire lagoon volume would recirculate through sediments every few days. These recirculation fluxes may therefore have important implications for the functioning of this coastal lagoon, since they may enhance the exchange of oxygen, solutes and particleassociated compounds between sediment and the overlying water column [START_REF] Anschutz | Tidal sands as biogeochemical reactors[END_REF][START_REF] Huettel | Transport and degradation of phytoplankton in permeable sediment[END_REF][START_REF] Huettel | Flow-induced uptake of particulate matter in permeable sediments[END_REF][START_REF] Jahnke | Organic matter remineralization and porewater exchange rates in permeable South Atlantic Bight continental shelf sediments[END_REF]. Accurately evaluating these recirculation fluxes is therefore required to understand the role that this process may play on the biogeochemical cycles of lagoon water and sediments.

Extension to Other Radionuclide Tracers

Within the last few years, a number of studies have used the 224 Ra/ 228 Th ratio in coastal sediments to calculate the 224 Ra flux and the corresponding water flux across the sediment-water interface [START_REF] Cai | Measurement of 224Ra:228Th disequilibrium in coastal sediments using a delayed coincidence counter[END_REF][START_REF] Cai | 224Ra:228Th disequilibrium in coastal sediments: Implications for solute transfer across the sediment-water interface[END_REF][START_REF] Cai | Using 224Ra/228Th disequilibrium to quantify benthic fluxes of dissolved inorganic carbon and nutrients into the Pearl River Estuary[END_REF]. Unlike 222 Rn, 224 Ra and 228 Th will partition onto the solid phase, and so the 224 Ra deficit must be calculated from the total exchangeable 224 Ra and 228 Th activities, rather than from 224 Ra and 228 Th dissolved in pore water. However, calculation of water fluxes from 224 Ra and 228 Th profiles also requires information on 224 Ra activities in the dissolved phases, as only this component will be transported with the water. In addition, the partitioning of 224 Ra between water and solid phases is largely dependent on the chemical composition of porewaters (e.g. salinity, pH, temperature, redox potential) and the characteristics of the sediments (e.g. grain size, fraction of exchangeable sites, content of iron and manganese), requiring an appropriate characterization of the 224 Ra solid-solution partitioning coefficients [START_REF] Beck | Controls on solid-solution partitioning of radium in saturated marine sands[END_REF]Gonnea et al., 2008). Finally, it should be noted that partitioning to the solid phase means that depletion of 224 Ra in pore water will be much shallower than for 222 Rn (as depleted 224 Ra in pore water will be replaced by 224 Ra released from the sorbed phase), and so this requires much finer resolution sampling.

Recent studies using the 224 Ra/ 228 Th method have been mainly in finer grained sediments than those that have used the 222 Rn method, and in these environments diffusion often forms a significant component of the tracer flux. The advection cycling model presented in this paper is less amenable to situations in which molecular diffusion is a significant component of the tracer flux. However, notwithstanding the above limitations, in advection-dominated systems, numerical models similar to that presented in this paper could potentially be used to estimate the variation in water flux with depth based on measured 224 Ra profiles.

CONCLUSIONS

We simulated radon porewater profiles using an advective cycling numerical

  c is concentration, Δc is change in concentration, 𝑞 ! and 𝑞 ! are the vertical and horizontal water fluxes, Δz is the vertical cell dimension, Δt is the temporal discretization and the superscripts 𝑎 and 𝑏 denote the fluxes and concentrations in the downward-flow and upward-flow cells respectively (symmetry of the recirculation cell requires that 𝑞 ! ! = -𝑞 ! ! ). After a period of time (tr/2), the flow reverses, and so the directions of all arrows shown in Figure 1 reverse. The calculations are then repeated, with superscripts 𝑎 and 𝑏 switched in Equations 6 and 7. This cycle is then repeated. Because the direction of flux changes during the recirculation cycle, the mean upward and downward water flux at each depth across a complete recirculation cycle (𝑞 ! ) is calculated by dividing 𝑞 ! by two (i.e. 𝑞 ! = 2𝑞 ! ).

  The mean upwelling or downwelling water flux across the sediment -water interface in Pz1 is 𝑞 ! 0 = 0.25 m d -1 (3 × 10 -6 m s -1 ). The radon flux from the surface into the underlying cell is 𝑞 ! 0 𝑐 0 , which is equal to 7.5 Bq m -2 d -1 . The radon concentration in this second cell (immediately below the surface water layer) is 240 Bq m -3 , and so there is also an upward radon flux into the surface water of 𝑞 ! 0 𝑐 1 = 61 Bq m -2 d -1 . The net radon flux into the surface water is therefore 𝑞 ! 0 𝑐 1 -𝑐 0 = 54 Bq m -2 d -1 . This is essentially the same as the radon flux for Pz1 calculated from Equation2(58 Bq m -2 d -1 ).

Figure 5

 5 Figure 5 depicts the sensitivity of the radon profiles to variations in surface water-

  Figure 1. Schematic representation of advective mixing cell model. Arrows denote

Figure 2 .

 2 Figure 2. Observed radon profiles at two different locations within La Palme

Figure 3 .

 3 Figure 3. (a) Mean vertical water velocity (upward or downward), as a function of

Figure 4 .

 4 Figure 4. Sensitivity of radon profiles to recirculation times (tr). Upwelling and

Figure 1 .

 1 Figure 1. Schematic representation of advective mixing cell model. Arrows denote

Figure 2 .

 2 Figure 2. Observed radon profiles at two different locations within La Palme

Figure 3 .

 3 Figure 3. (a) Mean vertical water velocity (upward or downward), as a function of

Figure 4 .

 4 Figure 4. Sensitivity of radon profiles to recirculation times (tr). Upwelling and

Figure 5 .

 5 Figure 5. Sensitivity of radon profiles in Pz1 to (a) surface water -sediment

Figure 6 :

 6 Figure 6: Simulated 𝑞 ! at different depths from 4 different radon profiles that

Table 1 .

 1 Measured values of porewater salinity and radon concentration. 867Uncertainty for radon concentration represents one standard deviation. 868

		Pz1			Pz2
	Depth	Salinity	222 Rn	Depth	Salinity	222 Rn
	(m)	(g/kg)	(Bq m -3 )	(m)	(g/kg)	(Bq m -3 )
	0	30.6	25 ± 10	0	37.2	89 ± 28
	0.05	32.9	933 ± 152	0.10	38.1	1238 ±
	0.10	30.8	1753 ± 218	0.15	38.1	1807 ±
	0.15	31.3	2895 ± 280	0.20	38.3	1953 ±
	0.20	30.8	3477 ± 327	0.30	38.6	2228 ±
	0.30	30.7	4258 ± 393	0.50	38.9	1925 ±
	0.50	37.7	4985 ± 496	0.80	51.6	2665 ±
	0.80	84.5	4022 ± 478	0.13	62.2	2600 ±

ACKNOWLEDGEMENTS

This research is a contribution to the ANR @RAction chair medLOC (ANR-14-ACHN-0007-01) and Labex OT-Med (ANR-11-LABEX-0061, part of the ''Investissements d'Avenir" program through the A*MIDEX project ANR-11-IDEX-0001-02), funded by the French National Research Agency (ANR). PC acknowledges support from IméRA (Institute of Advacned Studies), Aix-Marseille Université (Labex RFIEA and ANR "Investissements d'avenir"). This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 748896. We thank Jordi Garcia-Orellana and Joan Manuel Bruach from the Universitat Autònoma de Barcelona for radon sample analyses.

obtained for larger cycles. Note also that the inflection point in the downwelling profile is similar to the estimated water flux (qv(0)) multiplied by the period of downwelling (tr/2). Other model parameters are the same as in Figure 3.