C Hubert

K J Kubiak

M Bigerelle

A Dubois

L Dubar

Identification of lubrication regime on textured surfaces by multi-scale decomposition

Keywords: Roughness, Friction, Curvature radius, Fractal

This paper focuses on the analysis of surfaces resulting from thickness reduction of aluminium strips, provided with lubricant reservoirs, and deformed with di erent lubricant viscosities and drawing velocities, by a specific experimental drawing process. During reduction, the lubricant is pressurized and may escape from its initial cavities and supply the neighbouring ones. The nature of the lubrication regime is thus locally changed and may vary from boundary to hydrodynamic. The deformed surfaces are measured by means of a Vertical Scanning Interferometer then analysed in terms of arithmetic roughness and developed roughness profile length. Since these two parameters are found not accurate enough to identify regions where di erent lubrication phenomena occurred, a method based on roughness peaks and valleys curvature radii estimation, recently developed by some of the present authors, is applied.

Introduction

Metal forming operations generally require mixed lubrication regime to reduce the friction at the tool/workpiece interface and to improve the visual aspect and roughness of the final surfaces [START_REF] Bloor | An elasto-plasto-hydrodynamic lubrication analysis of the plane strain drawing process[END_REF]. The peculiarity of this lubrication regime lies in the fact that only the micro-cavities, inherent from surfaces roughness [START_REF] Kudo | Aninvestigationintoplasto-hydrodynamiclubrication with a cold sheet drawing test[END_REF], are filled with lubricant during the contact, and thus the peaks of both surfaces are in dry contact. In plane strain processing, when plastic deformation occurs, the cavities filled with lubricant are pressurized and release the lubricant towards the drawing direction or in the opposite direction, according to the involved phenomena: Micro Plasto HydroStatic Lubrication (MPHSL, governed by the pressure gradient between the cavity and the contact pressure) or Micro Plasto HydroDynamic Lubrication (MPHDL, mainly governed by the relative speed between the tool and the worked surface), respectively. The result of this lubricant escape is a local change of lubrication regime that may become hydrodynamic rather than the initially expected, mixed lubrication regime.

The concept of Micro Plasto Hydrodynamic Lubrication (MPDHL) was first proposed by Mizuno and Okamoto [START_REF] Mizuno | E ects of lubricant viscosity at pressure and sliding velocity on lubricating conditions in the compression-friction test on sheet metals[END_REF] and later by Kudo et al. [START_REF] Kudo | Aninvestigationintoplasto-hydrodynamiclubrication with a cold sheet drawing test[END_REF]. Azushima et al. [START_REF] Azushima | Direct observation of contact behaviour to interpret the pressure dependence of the coe cient of friction in sheet metal forming[END_REF] and Bech et al. [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF] proved it in plane strain drawing of aluminium strips on which pyramidal indentations were manufactured and filled with lubricant. The lubricant dynamic behaviour was observed through a transparent glass die and recorded by means of a video camera. The aim of the latter study was to determine the drawing parameters, such as the lubricant viscosity, drawing speed, back tension or friction coe cient, that lead to MPHSL and/or MPHDL. This study also showed strong changes in the worked surface roughness where MPHSL and MPHDL occur, but has not been investigated further. More recently, Dubar et al. [START_REF] Dubar | Analysis of fluid lubrication mechanisms in metal forming at mesoscopic scale[END_REF] carried out similar kind of experiments, with the same testing device as [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF], but with triangular sectioned grooves as lubricant reservoirs, machined by means of Electro Discharge Machining (EDM). The aim of this technique was to avoid metal banks due to the indentation process, to avoid any lubricant escape perturbation due to these banks and to confer true plane strain properties to the analysis, including the lubricant escape. The drawing parameters investigated in this study were the drawing speed, the lubricant viscosity and the pockets angle to the edge, with goal to model numerically the lubricant reservoirs depth reduction due to plastic deformation and subsequent fluid escape during drawing, by means of a Fluid Structure Interaction model. Other studies focusing on related phenomena in drawing process like friction, lubrication, roughness etc. have been carried out [START_REF] Lee | Modeling of the friction caused by lubrication and surface roughness in sheet metal forming[END_REF][START_REF] Gong | Studies on the friction and lubrication characteristics in the sheet steel drawing process[END_REF][START_REF] Christensen | Stochastic models for hydrodynamic lubrication of rough surfaces[END_REF][START_REF] Le | Measurements of friction in strip drawing under thin film lubrication[END_REF][START_REF] Costa | E ects of die surface patterning on lubrication in strip drawing[END_REF]. In the present study, the experimental results produced by Dubar et al. [START_REF] Dubar | Analysis of fluid lubrication mechanisms in metal forming at mesoscopic scale[END_REF] are used to analyse the e ects of lubricant escape on the final roughness of the worked surfaces.

As classical tribological parameters such as the arithmetic roughness (Ra) or the developed profile length (RL0) do not allow clear identification of regions where di erent phenomena occur [START_REF] Bigerelle | Roughness signature of tribological contact calculated by a new method of peaks curvature radius estimation on fractal surfaces[END_REF], the method based on roughness peaks curvature estimation, which plays an important role in contact pressure, recently developed by Bigerelle et al. [START_REF] Bigerelle | Roughness signature of tribological contact calculated by a new method of peaks curvature radius estimation on fractal surfaces[END_REF] is applied. This method considers the specimens surfaces as fractal surfaces and thus allows to perform a multiscale decomposition of this radius of curvature. Even if the calculation of the radius of curvature requires smooth surfaces, it is not always the case when dealing with fractal surfaces, as shown by Mandelbrot and Wallis [START_REF] Mandelbrot | Computer experiments with fractional gaussian noises. part 1, sample graphs, averages and variances[END_REF]. Moreover, as all metric parameters relative to fractal curve depend on the measurement scale, it becomes difficult to give a physical sense of the local radius of curvature for fractal surfaces. And if in tribology contact fractal surfaces are often used to avoid sensitivity of the measurement scale, this approach is used here to identify the regions where di erent contact conditions occurred, especially a local change in lubrication regime due to either MPHSL or MPHDL.

The present paper is divided into three parts. Firstly, the experimental tests carried out are described: the testing device functioning is detailed as well as the specimens preparation and properties, and the testing conditions. Then the methodology for the acquisition of the worked specimens topographical information is explained. Secondly the methodology to perform the whole multiscale decomposition of the radius of curvature along a given topographical profile is described. Then, based on the prior multiscale decomposition, the analysis of the results is carried out to quantify the e ects of the experimental process parameters (i.e. the lubricant viscosity, the drawing velocity and the pockets angle) on the radius of curvature measured on the final surfaces. Finally the link between this radius of curvature and the lubrication regime is presented, with comparisons to regions outside of the pockets row and to initial profiles. The results are then discussed and proposals for future works are given.

Experimental tests and surfaces topography acquisition

Testing device and design of experiments

The experimental testing device used in this study has been developed by Bech et al. [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF] for thickness reduction of aluminium strips. The strips are drawn by means of a hydraulic jack between two dies: the upper one is made of hardened glass so the lubricant behaviour can be observed during thickness reduction, and the lower die is made of steel with an inclined plane with angle =3

喰 to apply the thickness reduction. The hardened glass is used as delivered in shape of circular disc of dimensions 挨50 × 11 mm, and the lower die is 80 mm long and 50 mm wide, with polished contacting surface to reduce friction. A high speed video camera is used to record the lubricant behaviour for subsequent analysis, i.e. identification of the occurrence of MPHSL and MPHDL, location of the onset and direction of escape or front wave speed. The testing device is illustrated in figure 1.

Figure 1: Experimental testing device with its intrinsic parameters [START_REF] Kubiak | Dynamic evolution of interface roughness during friction and wear processes[END_REF].

The provided strips dimensions are l0 × L0 × t0 = 450 × 20 × 2 mm and their material is a semi-hardened aluminium AISI 1050 H24. Macro lubricant pockets are manufactured by 3 means of Electro Discharge Machining to reduce the local residual stresses and to avoid the formation of banks that would change the tribological contact properties and fluid escape, as observed in a prior study [START_REF] Hubert | Numericalsimulationoflubrication mechanisms at mesoscopic scale[END_REF]. The shape of the lubricant pockets is a groove normal to the drawing direction with triangular cross section. They are 10 mm long and 1 mm wide, with an angle to the edge =5 喰 for the first set of specimens, and = 10 喰 for the second one. The spacing between each cavity centre is 2 mm, so that the plateaus are 1 mm wide. The strips are obtained from cold rolling and their raw roughness is preserved all along the specimens preparation process. An illustration of the manufactured strips and an illustrative picture are given in figure 2. In the experiments carried out for this study several testing conditions are applied to analyse the sensitivity of the lubrication phenomena described above, i.e. MPHSL and MPHDL: the pockets angle to the edge applied during the EDM of the specimens ( =5 喰 and = 10 喰 ), the lubricant, which is a pure mineral oil with two viscosities ( = 60 cSt and = 660 cSt at 40 喰 C) and the drawing velocity (Vs =5 mms 梓1 and Vs = 50 mms 梓1 ). All these parameters are crossed together so that eight experimental conditions are tested, each testing condition being tripled resulting in twenty four tests. The design of experiments is summed up in table 1 along with the name assigned to each testing condition (e.g. 5DTC1 means 5 Degrees, Testing Condition 1).

During testing global parameters are recorded such as the drawing speed, drawing force, and the lubricant dynamic behaviour by means of a high speed video camera set to a frame rate of fr = 125 fps when drawing at Vs =5 mms 梓1 , and fr = 1500 fps at Vs = 50 mms 梓1 .

Deformed surfaces topography acquisition

The surfaces topography were acquired using a Zygo NewView 7300 [START_REF]Zygo NewView[END_REF] profilometer which is a Vertical Scanning Interferometer (VSI). The measurement principle is that unfiltered white light beam is split in two: half of the beam is directed through a microscope lens and reflected from the surface, and the other half is reflected from the reference mirror. When reflected beams combine together they produce interference fringes, where the best-contrast fringe occurs at best focus. In VSI mode the objective moves vertically to scan the surface at various heights. A 3D surface is then reconstructed by analysis of fringes at every pixel. VSI mode uses an algorithm to process fringe modulation data from the intensity signal to calculate the surface heights. The obtained resolution will therefore depend on the precision of the z-axis positioning. In the case of the Zygo NewView 7300 instrument a piezoelectric stage with range of 160µm is used to refine the height resolution going down to 0.1 nm. The spatial resolution depends on the camera size and the used lens; in our case the surfaces were acquired then stitched together to cover an area of 8 × 2 mm 2 with spatial resolution (x, y) of 1.1µm [START_REF] Kubiak | Dynamic evolution of interface roughness during friction and wear processes[END_REF]. The location of the measurements is chosen so that both deformation and lubrication regimes are stabilized. It starts at the 15th pocket, in the middle of the strip as illustrated in figure 3. been flattened and the resulting plateaus surface is smoother than the original surface. It can also be observed that the shape of the lubricant pockets has changed after processing, and two new regions can be identified, labelled 2 and 4 (figure 4b). Initially, those surfaces belonged to the lubricant pockets and, due to elastoplastic deformation, they now form part of a plateau. It is however assumed that because of the large amount of lubricant in this region, the contact was only partial and thus the topography of regions 2 and 4 is similar to that of the pockets, labelled 1. In this kind of deformation process, the surface roughness will change leaving a specific signature of tribological process that can be revealed with an analysis of the surface roughness. However, direct analysis of 3D surface morphology is computationally more expensive and less robust than 2D analysis in the case of pure anisotropic surfaces. Therefore, to analyse surface roughness, a specific technique of transition from 3D surface to 2D surface profiles analysis has been developed. Given a 3D surface, it consists in generation of series of 2D profiles in the y-direction of this surface. Once the profiles are extracted, di erent parameters can be calculated from individual profiles such as the classical arithmetic roughness or the developed length of the roughness profile, as described in [START_REF] Bigerelle | Roughness signature of tribological contact calculated by a new method of peaks curvature radius estimation on fractal surfaces[END_REF] and illustrated in figure 5. Though, as described in [START_REF] Bigerelle | Roughness signature of tribological contact calculated by a new method of peaks curvature radius estimation on fractal surfaces[END_REF], using a simple amplitude parameter such as the arithmetic roughness (Ra) can not clearly distinguish all contacting regions. The developed length of roughness profile (RL0) is able to distinguish them with slight variations, but is not clear enough to perform a reliable analysis. Thus, an algorithm has been developed to estimate the radius of curvature of a fractal curve and to establish the dependency of this radius on the scale at which the surface is observed and its relationship with the fractal dimension of this surface. This work is fully detailed in [START_REF] Bigerelle | Roughness signature of tribological contact calculated by a new method of peaks curvature radius estimation on fractal surfaces[END_REF] and the main lines are given in section 3.1.

Multiscale decomposition of profile curvature radius

Methodology

According to Nowicki [START_REF] Nowicki | Multiparameter representation of surface roughness[END_REF], it is possible to study the roughness of a given surface by introducing a parameter called radius of asperity, which is determined by the relation rc =lx 2 Equation 1 is obtained by assuming that ly is su ciently small so that condition ly << lx is satisfied.

For fractal curves, rc depends on the scale at which the observation is made, and as we postulate that the curvature radius can be defined at a given scale, the parameter lx is assumed to be the measurement scale, calculated as follows: 1. A horizontal straight line is chosen at a given level h that crosses the measured profile and a set of values of lx is built with the intersections of the straight lines intercepting the profile. 2. For each value of lx the local maximum peak (i.e. the maximum value of the profile) is obtained, giving the corresponding value of ly.

3. The curvature radius of the peak rc is calculated according to equation 1. The multiscale decomposition will be performed by means of this scale law.

Application on an initial profile

The algorithm is introduced with a 2D profile while the acquired data are surfaces, and thus 3D. Meanwhile the experimental tribotest used for this study performs plane strain deformations, which means that the mechanical behaviour of the specimens is (theoretically) identical whatever the observation position in the width direction. As done for parameters Ra and RL0 (figure 4), the algorithm will be applied on all profiles in the width direction and theobtainedvalueofthecurvatureradiusofpeakswillbereportedona log rc(lx)= f(xprofile) graph for several values of the measurement scale log lx, ranging from 0.35 to 3.0. As discussed earlier, the radii of curvature of the profile peaks are an important indicator of contact severity and thus contact pressure. In the same way as developed, another important indicator may be the radius of curvature of valleys: it is clear that for soft to medium contact the curvature radius of peaks will change, but the curvature radius of valleys won't change significantly. On the contrary for severe contact conditions, the peaks will flatten leading to increased curvature radius, as well as the valleys. To compute the radius of the profile valleys, the same algorithm as for the peaks is applied after initial transformation of the surface z(x, y) into -z(x, y). Finally, still based on the radii of curvature of peaks and valleys, a third indicator will be useful for the interpretation of the results: the dimensionless ratio R between the radius of curvature of peaks and the radius of curvature of valleys, expressed as

R = r cP , (5) rcV
with rcP and rcV the radii of curvature of peaks and valleys, respectively. In logarithmic scale, as it will be used in the remainder of this paper, this indicator becomes

R 念 = log R = log rcP -log rcV . ( 6
)
Three main values can be considered for this indicator:

1. R 念 > 0: the radius of peaks is greater than that of valleys, meaning that the peaks have been eroded but not the valleys. The other way may also be considered, i.e. rcV becomes smaller without change for rcP , but it seems to be inconsistent since any process starts by eroding the peaks then the valleys, except if the valleys are filled with third body particles, created by the tool/workpiece friction and wear.

2. R 念 =0: means that both radii of curvature (peaks and valleys) are equal.

3. R 念 < 0: the peaks curvature radius has been reduced or that of the valleys has been increased.

These three indicators, namely the peaks curvature radius, the valleys curvature radius and the ratio between these last two curvature radii have been computed for the initial strip illustrated in figure 4a, and reported in figure 7, along with a top view of the surface topography and a 2D profile of this surface, taken in its middle. It can be seen on the graph of figure 7 that the curvature radii of peaks rcP and valleys rcV are very close, if not equal. The analysis and explanations of the di erent parameters is given below:

• Logically, regions that have been cold rolled have high curvature radii since cold rolling apply large plastic deformation to the strips leading to important flattening of both peaks and valleys. Moreover, the work rolls being much harder than the aluminium strip, they transfer their roughness, which is typical of rolls grinding, to the worked strips.

•

The regions a ected by the Electro Discharge Machining also have similar value of curvature radii of peaks and valleys, but this value is smaller than that encountered for the strips plateaus. The EDM actually removes aluminium particles almost uniformly instead of flattening the strip surface.

•

For both regions, the ratio R 念 of peaks and valleys curvature radii (equation 6) is thus close to zero, ranging from -0.2 to +0.2 for all measurement scales.

These values of the peaks and valleys curvature radii and ratio for this initial profile will be kept as reference values for the next analysis of the measurements performed after drawing.

Determination of the most relevant analysis scale

The major interest and key for the following study is to find the scale i.e. the lx value at which the ratio R 念 of peaks and valleys curvature radii must be analysed. To determine this scale, appropriate statistical analyses are performed. For all profiles, all values of individual peaks and valleys curvature radii are taken into account, giving 452, 861, 509 peaks and 452, 861, 510 valleys, so about 9 × 10 8 asperities. From these analyses, it is decided to test if the di erence between peaks and valleys is more significant for all measurement scale values lx and this for all maps and all x-positions. Then analyses of variance are proceeded to find where is the more relevant measurement scale lx from which di erence between peaks and valleys curvature are more significant and thus all cases (horizontal x-position and each sample). Only the most significant values are retained. The mean value of F (Fisher variate of the analysis of variance) variation versus the measurement scale lx (in logarithmic values) are plotted on the graph of figure 8. Higher F correspond to higher di erence between peak and valley curvatures. As it can be observed in figure 8, two maximum values are reached: the first one around log lx =0.35 (lx =2.24µm) and the second one around log lx =2.5 (lx = 315.23µm). These two measurement scale values are linked to the size of the asperities; the others represent the macrotopography linked to the structuration of the specimens. As a consequence, to analyse e ect of the topography, we will retain the measurement scale log lx =0.35.

Analysis of the deformed surfaces and identification of lubrication regimes

This section concerns the analysis of the specimens after drawing, for all configurations detailed in table 1. It is divided into 3 sub-sections: the first one is dedicated to the analysis of the e ect of the lubricant viscosity, the second one will concern the drawing velocity and the third sub-section will finally be aimed at the identification of the different lubrication regimes that are encountered in this specific kind of tribological contact, by comparisons with measurements performed out of the pockets rows. The first two sub-sections will also deal with the e ect of the pockets angle to the edge, which effects on the curvature radii are linked with that of the lubricant viscosity and the drawing velocity.

E ect of the lubricant viscosity

Generally, in mixed lubrication regime, increasing the lubricant viscosity allows to reduce the friction between the contacting surfaces and may favour a hydrodynamic lubrication regime.

The expected e ect of the high viscosity lubricant ( = 660 cSt) is thus smaller curvature radii for both peaks and valleys belonging to the strips plateaus (closer to the initial profiles) than for a low viscosity lubricant ( = 60 cSt). A comparison for testing conditions 5DTC1 and 5DTC2 is given in figure 9a and 9b, which correspond to low velocity tests with small and high lubricant viscosity, respectively. On this figure, it can be seen that, even if not constant, the value of log rcP is higher for the low viscosity lubricant (5DTC1, figure9a)andrangesin [1. Still focusing on the strips plateaus, the same observations can be made for testing 14 conditions 5DTC3 and 5DTC4, that involve a high drawing velocity and low and high lubricant viscosities, respectively. For these configurations the radius of curvature of peaks is log rcP ≈ 1.3 (rcP ≈ 19.95µm) for the low viscosity lubricant while log rcP ≈ 0.8 (rcP ≈ 6.31µm) when using a high viscosity lubricant. For the valleys, log rcV ranges in the interval [0.8; 1] (rcV 樺 [6.31, 10.00]µm)fortestswithlowviscositywhile log rcV ≈ 0.6 (rcV ≈ 3.98µm) for tests with high viscosity lubricant. Increasing the lubricant viscosity logically protects the contacting surfaces where the peaks contact pressure should be found smaller than for a low viscosity lubricant. All tests performed with pockets angle = 10 喰 (i.e. tests conditions names 10DTCx, x being the test combination) show the same trend: compared to the initial strips state (figure 7), the curvature radius of both peaks and valleys have much increased when using a low viscosity lubricant than when using a high viscosity lubricant. Meanwhile for = 10

喰 , theobservedcurvatureradiiareslightlyhigherthanforstripswithpocketsangle =5
喰 , and especially for the peaks curvature radius. The multiscale decompositions for two identical running conditions but with two di erent angles are given in figure 10a and10b, for =5 喰 and = 10 喰 respectively. plays an important role in the hydrodynamic part of contact lubrication. In their study, which is experimentally similar to that of Bech et al. [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF] and to the present one, specimens provided with pyramidal indentations and several angles to the edge were drawn. The determinations of the onset of lubricant escape from the pockets by MPHDL showed that it occurred earlier for small pyramid angles than for high angles, which ranges from =3.6 喰 to = 15 喰 .

Itthusleadstolesslubricantonthestripsplateausandmoreseveretool/workpiece contact.

The same phenomenon is observed in the present study and is highlighted by the curvature radii of peaks and valleys on the strips plateaus: MPHDL occurs later for pockets with = 10 喰 than for pockets with =5 喰 , with consequentially less lubricant on the strips plateaus, leading to more severe contact and thus, high curvature radii. Moreover, the value of log rcP is not only slightly higher for 10DTC2 than for 5DTC2, its profile is also di erent: it has a negative slope towards the +x-direction, that corresponds to the direction of MPHDL, while for the small pockets angle the profile of log rcP is almost flat. As illustrated in figure 11, it means that the left border of the strips plateaus have undergone severe contact due to late lubrication. MPHSL, which is less sensitive to pockets angle do occurred and had lubricated the right part of the strips plateaus. In the case of configuration 5DTC2, MPHDL occurred earlier and both MPHSL and MPHDL melted to result in su cient lubrication and almost uniform peaks radius of curvature on the strips plateaus.

E ect of the drawing velocity

The drawing velocity also plays an important role in contact lubrication: as stated by Sørensen et al. [START_REF] Sørensen | A basic study of the influence of surface topographyon mechanismsofliquid lubricationinmetal forming[END_REF] or Bech et al. [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF], it increases the pressure in the lubricant pocket which overcomes the contact pressure (also called sealing pressure) and thus helps the lubricant to be dragged out of the cavity. In the case of the experiments carried out in this study, this e ect can also be observed and, since a hydrodynamic lubrication regime tends to preserve the surfaces roughness, the observed final curvature radii are closer to the initial ones than for the tests performed at low velocity.

For strips with pockets angle =5 Again, for strips with pockets angle to the edge = 10 喰 , the same trend is observed but the e ect of the pockets angle on the MPHDL onset discussed in the previous section is reinforced. IftheincreaseofthedrawingvelocitytendstoincreasetheMPHDLphenomenon [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF], it also tends to reduce the extent of the MPHSL [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF]. The so-called lubricant forward escape [START_REF] Dubar | Analysis of fluid lubrication mechanisms in metal forming at mesoscopic scale[END_REF] will be less far towards the contact exit (-x-direction) than when drawing at low velocity, and thus the ejected lubricant will remain close to its initial pocket and only lubricateitsleftedge. Analogouslytotheillustrationoffigure11b,theresultisadescendant slope for the profile of the peaks radius of curvature according the +x-position on the strips plateaus. Figure 12 shows this phenomenon through the multiscale decompositions related to testing conditions 10DTC2 and 10DTC4. decreased (from log rcP ≈ 1.6, rcP ≈ 39.81µm at low velocity to log rcP ≈ 1.5, rcP ≈ 31.62µm), but the minimal value on the strips plateaus decreased with the increase of the velocity (from log rcP ≈ 1.1, rcP ≈ 12.59µm at low velocity to log rcP ≈ 0.8). Along with the decreased values of log rcV , the curvature radii measured for this testing condition are very close to the curvature radii of the initial profiles, especially in the right region of the strips plateaus. It shows that initial surfaces have not undergone severe contact conditions and thus contact pressure, contrarily to the opposite side of the plateaus.

Identification of lubrication regimes by means of the curvature radii

From a global point of view, the lubrication regimes involved in a strip drawing process can be estimated by means of process parameters and global measurements, such as the lubricant viscosity and the drawing forces, and their influence on the resulting friction coefficient. This influence has been described by Stribeck, and is commonly known as Stribeck's curve. In the case of this study, it is possible to plot such a curve, after friction coefficient calculation, and to estimate the global lubrication regime involved during the drawing operation of each specimen.

Two friction coefficients should ideally be identified, one for each tool surface. The testing device allows, however, only determination of one value based on measurements of the drawing force and inverse analysis using the slab method, implying that an average friction coefficient will be assumed.

Each testing configuration is plotted on the graph of figure 13. It must be noted that configuration 10DTC1 is not plotted since no stable regime (in terms of drawing force) has been reached during the tests, because of important sticking of aluminium strips to the glass die. Except for testing condition 5DTC1, it can be seen that the tests are quite reproducible, and all three repetitions for each condition present very small standard deviation (see error bars on the graph in Fig. 13). In the case of 5DTC1, it can be assumed that this condition results in boundary lubrication and increased coefficient of friction due to direct material/die contact.

Conditions 5DTC2, 5DTC3,10DTC2 and 10DTC3 are in a region between boundary and mixed lubrication regimes, and higher lubricant viscosity and increased drawing velocity in tests 5DTC4 and 10DTC4 result in fully mixed lubrication regime.

Taking into account the involved drawing velocities and the nature of lubrication during the drawing process, it is difficult to obtain a hydrodynamic lubrication regime and, as it can be noted in Fig. 13, the lowest values of friction coefficients represent a mixed lubrication regime.

As explained in introduction of this paper and as shown by Bech et al. [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF], local changes in lubrication regimes may occur because of MPHDL and MPHSL occurrence. Sections 4.1 and 4.2 highlighted these local changes in terms of peaks and valleys curvature radii, but these analyses are relative, based on comparisons of di erent testing conditions. They are important for the understanding of the contact lubrication and the effects of the drawing conditions on the lubrication, but the local lubrication regimes have not been estimated and this cannot be only based on the radius of curvature of peaks or valleys. Still based on comparisons, another region of the specimens can be useful to make a step forward in the analyses: regions outside of the pockets row, as illustrated in figure 14. This measurement location is not a ected by the pockets network and thus the measured surfaces have been deformed with constant lubrication regime. The measured surfaces are also longer than the plateaus, with the same width, that will allow to average the values of the peaks and valleys curvature radii. The averaged values of log rcP and log rcV and the value of the ratio of curvature radii R 念 are given in table 2, for all lubricant viscosities and drawing velocities combinations. Each combination is given a name similar to that used in table1, withouttheanglevaluesincethepocketsdonotplayanyroleinthesemeasurements.

A first look at the values of peaks curvature radius shows that they are consistent with the analyses of strips with pockets angle =5 喰 detailed in section 4.1. Concerning the curvature radius of valleys, they are slightly higher in regions out of the pockets row. This means that contact is more severe when no additional lubricant is provided in between the contacting surfaces. On the contrary, in the case of strips with pockets angle = 10 喰 , the curvature radii measured on the plateaus are higher than outside of the pockets row, which is supposed to be due to the sti pocket angle, that reduces the MPHDL e ects as discussed in sections 4.1 and 4.2.

As introduced in section 3.2, the ratio of peaks and valleys curvature radii in logarithmic scale (equation ( 6)) is here of great interest: it gives an indication on the morphology of the surface peaks and valleys that can be linked to the contact severity and thus, an indication of the lubrication regime can be given. As recalled in introduction, a mixed lubrication regime corresponds to few dry contacts between the tool/workpiece surfaces and a thin lubricant film in between these surfaces. Hydrodynamic lubrication is activated when there is almost no contact between the surfaces and, the opposite lubrication regime is the boundary regime, which takes place if the lubricant film is nearly broken. An illustration of these three main lubrication regimes is given in figure 15. 9,R 念 ranges between 0.45 and 0.65, with lower log rcV than for regions outside the pockets row. This means that the contact was less severe for this testing condition and confirms the protective e ect of the lubricant reservoirs. This contact lubrication would range between mixed and hydrodynamic (undeformed valleys and eroded peaks). More generally, all tested configurations with =5 喰 give the same range 2, both regions of the plateaus can be assigned a mixed lubrication regime while the global regime tends to hydrodynamic (Fig. 13 and Table 2). The combination of parameters R0 and log rcV allows to highlight a change in surface topography, and thus a change in local contact conditions due to local lubricant escape phenomena, namely MPHSL and MPHDL.

Conclusion

In this paper, a fractal approach of the curvature radius has been applied, on the basis of experimental tests and measurements. The aim was to analyse the effects of aluminium strips drawing process conditions on the final surfaces state and local changes in lubrication regime. After multiscale decomposition of the profile peaks and valleys curvature radii, the results corresponding to each testing configuration are compared to quantify the effects of change of lubricant viscosity, drawing velocity and lubricant pockets angle on the curvature radii of both peaks and valleys. These two parameters, in addition to a third one which is the ratio between the first two parameters, are used as indicator of contact severity: the higher the radius of curvature, the higher the contact pressure was during plastic deformation. Then, by comparison with regions of the specimens that have been deformed without any effect of the lubricant pockets, the different lubrication regimes encountered in the tribological contact have been determined qualitatively.

The approach applied in this paper seems to be quite relevant and accurate to distinguish the regions where different contact conditions occurred. Nevertheless, the authors have profoundly analysed each initial configuration, their profiles, roughness and curvature radii to eventually be more accurate in the identification of contact properties and lubrication regimes. It is also necessary to perform more basic experiments to establish more solid grounds for further investigations.

Figure 2 :

 2 Figure 2: Groove shaped pockets pattern (a) and picture (b) of two EDM processed strips.

Figure 3 :

 3 Figure 3: Illustration of measurements location and extent.

Figure 4 :

 4 Figure 4: Illustration of strips surfaces before (a) and after (b) drawing.

Figure 5 :

 5 Figure 5: Roughness parameters Ra and RL0 calculated on individual profiles, extracted in y-direction from the 3D surface [12].

  /8ly (1) where rc is the radius of a circle with centre O, passing by a crest A and by two points B and C of the profile. lx corresponds to the distance between points B and C, ly is the distance between line (B, C) and point A, andthedistance rc -ly correspondstoth edistancebetwe en line (B, C) and the circle centre O. ly is considered su ciently small so that segment [A, O] is supposed to be perpendicular to segment [B, C], in its middle I. These parameters are illustrated in figure 6 on a regular profile.

Figure 6 :

 6 Figure 6: Definition of parameters used to calculate the local radius of curvature (after [12]).

4 .

 4 The sequence[START_REF] Bloor | An elasto-plasto-hydrodynamic lubrication analysis of the plane strain drawing process[END_REF][START_REF] Kudo | Aninvestigationintoplasto-hydrodynamiclubrication with a cold sheet drawing test[END_REF][START_REF] Mizuno | E ects of lubricant viscosity at pressure and sliding velocity on lubricating conditions in the compression-friction test on sheet metals[END_REF] is repeated m times for m other horizontal straight lines chosen randomly.Experimentally, the fractal dimension is obtained as a slope by fitting the discretised data of the function log rc(lx)= f(log lx) obtained by the algorithm described above. If the regression line fit well (usually with a coefficient of correlation R ranging in [0.95, 1]), the experimental data then allow writing the scale law rc(lx)= lx .[START_REF] Azushima | Direct observation of contact behaviour to interpret the pressure dependence of the coe cient of friction in sheet metal forming[END_REF]

Figure 7 :

 7 Figure 7: Multiscale decomposition of peaks and valleys curvature radii and parameter R 念 for an initial strip with pockets angle of 10 丗 .

Figure 8 :

 8 Figure 8: Determination of the most relevant scale to observe the multiscale decomposition.

Figure 9 :

 9 Figure 9: Multiscale decomposition of peaks and valleys curvature radii for conditions 5DTC1 (a) and 5DTC2 (b).

Figure 10 :

 10 Figure 10: Multiscale decomposition of peaks and valleys curvature radii for conditions 5DTC2 (a) and 10DTC2 (b).

Figure 11 :

 11 Figure 11: Illustration of the e ect of the pockets angle to the edge on the evolution of the peaks curvature radius of the strips plateaus on configurations 5DTC2 (a) and 10DTC2 (b).

Figure 12 :

 12 Figure 12: Multiscale decomposition of peaks and valleys curvature radii for conditions 10DTC2 (a) and 10DTC4 (b).

Figure 13 :

 13 Figure 13: Evolution of the global friction coe cient according to the drawing velocity and the lubricant viscosity, for each test configuration.

Figure 14 :

 14 Figure 14: Location of measurements performed outside of the pockets row.

Figure 15 :

 15 Figure 15: Illustration of the main three lubrication regimes.

:

  close to 0.5 for the low viscosity lubricant, smaller values for the high viscosity lubricant, around R 念 =0.35 with values of log rcV close to the initial ones. It is more interesting to observe the strips plateaus for testing configurations with = 10 喰 , for which the profiles of curvature radii were not uniform. Taking the example of testing condition 10DTC2 (illustrated in figure 12a), on the strips plateaus R 念 varies in the range [0.32, 0.42], while for 10DTC4 (illustrated in figure 12b) R 念 varies in the range [0.25, 0.45] with smaller values of log rcV . The postulated reason for this change in curvature radii was the e ect of the sti pockets angle limiting the hydrodynamic lubrication, leading to severe contact on the left side of the plateaus (high values of R 念 with high values of log rcV ), and soft, well lubricated contact on the right side of the plateaus (figure 11b, small values of R 念 with values of log rcV close to the initial ones). In terms of ratio of curvature radii and for condition 10DTC4, the left regions of the plateaus have a ratio R 念 ≈ 0.45 with high values of log rcV , while the opposite side, where MPHSL has lubricated the contact, have a ratio R « 0.25. According to the values of Table

Table 1 :

 1 Testing conditions and corresponding condition name.

  喰 , log rcP ≈ 1.3 (rcP ≈ 19.95µm, while 1.3 ≤ log rcP ≤ 1.8, 19.95 ≤ rcP ≤ 63.1µm at low velocity) for the low viscosity lubricant and is reduced to log rcP ≈ 0.8, rcP ≈ 6.31µm(log rcP ≈ 1.1, rcP ≈ 12.59µm at low velocity) with the high viscosity lubricant, which is consistent with the analysis given in the previous section. The same trend is observed for the curvature radius of valleys. log rcV now ranges in [0.8, 1] ([6.31, 10.00]µm) instead of [0.8, 1.3] ([6.31, 19.95]µm) at low velocity and low viscosity, and it is reduced to log rcV ≈ 0.6 (rcV ≈ 3.98µm) with the high viscosity lubricant. This last value is very close to that of the initial profile, illustrated in figure7, and clearly highlights the protective feature of hydrodynamic lubrication.

Table 2 :

 2 Averaged values of peaks and valleys curvature radii for all lubricant viscosities and drawing velocities combinations.

This is an pre-print of the paper published by Elsevier in Tribology International, Vol. 82 part B (2015) p. 375-386 kris@kubiak.co.uk https:/ / doi.org/ 10.1016/ j.triboint.2014.04.035

Acknowledgements

Authors would like to thank Professor Niels Bay from Technical University of Denmark for providing experimental support for tribological drawing tests and stimulating scientific discussions. ThepresentresearchworkhasbeensupportedbytheInternationalCampusonSafetyand Intermodality in Transportation, the Nord-Pas-de-Calais Region, the European Community, the Regional Delegation for Research and Technology, the Ministry of Higher Education and Research, and the National Centre for Scientific Research. The authors gratefully acknowledge the support of these institutions.


Sudarat Chadsuthi
email: sudaratc@nu.ac.th

Karine Chalvet-Monfray

Suchada Geawduanglek

Phrutsamon Wongnak

& Julien Cappelle

Spatial-temporal patterns and risk factors for human leptospirosis in Thailand, 2012-2018

Leptospirosis is a globally important zoonotic disease. The disease is particularly important in tropical and subtropical countries. Infections in humans can be caused by exposure to infected animals or contaminated soil or water, which are suitable for Leptospira. To explore the cluster area, the Global Moran's I index was calculated for incidences per 100,000 population at the province level during 2012-2018, using the monthly and annual data. The high-risk and low-risk provinces were identified using the local indicators of spatial association (LISA). The risk factors for leptospirosis were evaluated using a generalized linear mixed model (GLMM) with zero-inflation. We also added spatial and temporal correlation terms to take into account the spatial and temporal structures. The Global Moran's I index showed significant positive values. It did not demonstrate a random distribution throughout the period of study. The high-risk provinces were almost all in the lower north-east and south parts of Thailand. For yearly reported cases, the significant risk factors from the final bestfitted model were population density, elevation, and primary rice crop arable areas. Interestingly, our study showed that leptospirosis cases were associated with large areas of rice production but were less prevalent in areas of high rice productivity. For monthly reported cases, the model using temperature range was found to be a better fit than using percentage of flooded area. The significant risk factors from the model using temperature range were temporal correlation, average soil moisture, normalized difference vegetation index, and temperature range. Temperature range, which has strongly negative correlation to percentage of flooded area was a significant risk factor for monthly data. Flood exposure controls should be used to reduce the risk of leptospirosis infection. These results could be used to develop a leptospirosis warning system to support public health organizations in Thailand.

 .

 .

 . A study in Nan province of Thailand found that the open habitat near rivers or in rice fields that were prone to OPEN

flooding, were related to human leptospirosis infections [START_REF] Della Rossa | Environmental factors and public health policy associated with human and rodent infection by leptospirosis: a land cover-based study in Nan province, Thailand[END_REF] . The activities of rice farming are also associated with an increase of leptospirosis infections in other countries, for example, in Iran [START_REF] Sahneh | Investigation of risk factors associated with leptospirosis in the North of Iran (2011-2017)[END_REF] or Tanzania [START_REF] Maze | Risk factors for human acute leptospirosis in northern Tanzania[END_REF] . Thus far, to our knowledge, there have been few studies on the impact of rice production activities on leptospirosis infections.

Spatial-temporal variation of human leptospirosis incidences have also been investigated to identify high-risk areas, where preventive measure should be implemented [START_REF] Luenam | Modelling and analyzing spatial clusters of leptospirosis based on satellite-generated measurements of environmental factors in Thailand during 2013-2015[END_REF][START_REF] Baquero | Spatiotemporal dynamics and risk factors for human Leptospirosis in Brazil[END_REF][START_REF] Dhewantara | Geographical and temporal distribution of the residual clusters of human leptospirosis in China, 2005-2016[END_REF][START_REF] Nunez-Gonzalez | Spatial patterns of leptospirosis in Ecuador, 2013-2018[END_REF] . The study of leptospirosis spatial-temporal patterns may help to prioritize the prevention and surveillance of leptospirosis in high-risk areas. Also, risk factors should be identified in order to control human leptospirosis.

In this work, we established a large-scale dataset of human leptospirosis in Thailand over seven years from 2012 to 2018. The objectives of this study are (1) to investigate the spatial-temporal patterns of leptospirosis incidence, (2) to determine the high-risk and low-risk provinces across Thailand, and (3) to examine the effect of rice-production-related factors, remotely sensed environmental factors, and climatic factors on leptospirosis risk in humans using a generalized linear mixed model (GLMM).

Materials and Methods

Data collections. Human leptospirosis case reports were retrieved from the database of national disease surveillance (report 506), Bureau of Epidemiology, Department of Disease Control, Ministry of Public Health, Thailand 19 . Most positive cases were suspected leptospirosis cases, based on the clinical diagnosis made by attending physicians and were mainly reported from public hospitals. In this research, we analysed all reported cases from 2012 to 2018 at the province level.

The percentage of flooded area was obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) of the Terra satellite (Surface Reflectance 8-Day L3 Global 500 m SIN Grid V005 (MOD09A1)). Using the surface spectral reflectance, band 4 (green) and band 7 (infrared), which measured at ground level in the absence of atmospheric scattering or absorption, the modified normalized difference water index (MNDWI) was calculated [START_REF] Ledien | Assessing the performance of remotely-sensed flooding indicators and their potential contribution to early warning for leptospirosis in Cambodia[END_REF][START_REF] Xu | Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery[END_REF] . To calculate the percentage of flooded area, the number of pixels for which the MNDWI value was greater than or equal to zero, in each province, for each month and for each year, were counted.

To reflect the vegetation intensity, the normalized difference vegetation index (NDVI) was calculated from the value of red band (wavelength: 620-670 nm) and near-infrared (NIR) band (841-876 nm) [START_REF] Hashim | Urban vegetation classification with NDVI thresold value method with very high resolution (VHR) PLEIADES Imagery[END_REF] . The vegetation area was classified into two types: low vegetation (such as shrub and grassland) and high vegetation (such as temperate and tropical urban forest). The vegetation type was calculated from the percentage of NDVI pixels that had a value within 0.2-0.5 and 0.5-1, respectively [START_REF] Hashim | Urban vegetation classification with NDVI thresold value method with very high resolution (VHR) PLEIADES Imagery[END_REF] .

The monthly values of rainfall, minimum and maximum temperature, and soil moisture were obtained from TerraClimate [START_REF] Abatzoglou | a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015[END_REF] , calculated over a grid of approximately 4 km 2 , and aggregated at the province level. The annual average values were calculated from monthly data. We also added the temperature range (i.e., the average difference between maximum and minimum temperature) to test the model. Due to the relation of temperature range with cloud cover [START_REF] Yong | Variability of diurnal temperature range in east Africa during 1921-2010[END_REF] , we hypothesize that this parameter may be used as a proxy for rainfall and flooding.

Elevation data was derived from the NASA Shuttle Radar Topographic Mission (SRTM) 90 m Digital Elevation Data [START_REF]SRTM 90m DEM Digital Elevation Database[END_REF] as it may drive leptospirosis infection [START_REF] Dhewantara | Geographical and temporal distribution of the residual clusters of human leptospirosis in China, 2005-2016[END_REF] . The average elevation at the province level was used in this analysis. Slope data was calculated using QGIS version 3.6.0 from the elevation data in units of degree.

The human population data was obtained from the WorldPop database (http:// www. world pop. org) and summarised as human population density per year for each province. Livestock population density (buffalo, cattle, and pigs) were obtained from the Information and Communication Technology Center (ICT), Department of Livestock Development of Thailand at the district level [START_REF]of Thailand. Number of farmer households and livestock classified by species in Thailand[END_REF] .

In this work, we also incorporated rice production data such as rice crop arable area, rice cultivated area and rice yield, which may associate with leptospirosis occurrences. In Thailand, most cases were identified in rice farmers [START_REF] Hinjoy | Epidemiology of leptospirosis from Thai national disease surveillance system[END_REF] . Rice farmers typically plant rice in two seasons, i.e., the in-season rice is called the primary rice crop and the off-season rice is called the secondary rice crop. The primary and secondary rice season data was obtained from the Office of Agricultural Economics, Thailand (http:// www. oae. go. th/). The descriptions of the variables are presented in Table 1.

Spatial-temporal autocorrelation.

To determine the spatial autocorrelation, the annual incidence cases and monthly incidence cases per 100,000 population were used to calculate the Global Moran's I index [START_REF] Moran | The interpretation of statistical maps[END_REF] at the province level. The Moran's I value, which ranges from -1 to + 1, is used to classify the spatial clustering. A Moran's I value of more than zero indicates positive autocorrelation, and a value of less than zero indicates negative autocorrelation. When the Moran's I value is close to zero, it indicates a random distribution. The value is calculated from the spatial weight matrix, which is assigned from the nearest neighbour province. To test the significant level, we used a Monte-Carlo simulation method with 999 permutations. A significance of less than 0.05 is used to consider the cluster or no autocorrelation. We used R package 'spdep' to calculate the Moran's I value [START_REF] Bivand | Spatial dependence: Weighting schemes, statistics[END_REF] .

We also determined the local indicators of spatial association (LISA) at the province level. LISA was used to describe leptospirosis clusters with hot spots (High-High), cold spots (Low-Low) and spatial outliers (High-Low and Low-High) [START_REF] Anselin | Local indicators of spatial association-LISA[END_REF] . Hot spots or High-High provinces indicate provinces that had a high leptospirosis incidence rate and were surrounded by provinces with high incidence rates. Cold spots represent provinces with low incidences surrounded by other low incidence provinces. Spatial outliers represent high incidence provinces surrounded by low incidence provinces (High-Low) and low incidence provinces surrounded by high incidence provinces (Low-High). LISA was calculated using GeoDA version 1.18 software (https:// geoda center. github. io/) [START_REF] Anselin | An introduction to spatial data analysis[END_REF] . The spatial weights file was created based on an inverse distance function with a distance band of 150 km. We considered a p-value less than 0.05 as a statistically significant result using 999 permutations.

Statistical model.

To study the spatial and temporal distribution of leptospirosis, we used a generalized linear mixed model (GLMM). Mixed effect models were used to account for the spatial dependence using province ID and for the temporal dependence, using time as random effect variables. The response variable of this model was the number of reported cases. Due to the data containing many zeros, the model would underpredict the cases. In this work, we used a zero-inflated negative binomial distribution (ZINB), which is an extension of the Poisson regression and captures the excess number of zeros [START_REF] Brooks | glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling[END_REF][START_REF] Brooks | Modeling zero-inflated count data with glmmTMB[END_REF] . To fit ZINB, we used the R package 'glm-mTMB' [START_REF] Magnusson | glmmTMB: Generalized linear mixed models using template model builder[END_REF] . We also added spatial and temporal correlation terms to consider the spatial and temporal structures (Table 1).

All variables were scaled by subtracting their mean and dividing by the standard deviation before investigating the effect of factors. We first selected the variables by studying the cross-correlation. For pairs of variables that have highly correlated, (more than |0.5|), only one variable was selected based on the minimum AIC value of the univariate model. Then, multivariate models, including all selected variables, were examined using the dredge function from the R package 'MuMIn' [START_REF] Barton | Multi-model inference[END_REF] . This function generates the 2 N model table with subsets of fixed effect terms, where N is the number of variables. The best models are sorted based on the criterion (AIC). The final model was the model with the least number of variables in the set of models with a ΔAIC < 2 in comparison with the model of the lowest AIC.

Ethics statement. All information of leptospirosis surveillance data were collected from the Thai Ministry

of Public Health. This study was approved by the Institutional Review Board of Naresuan University (P10003/64). The need of informed consent was waived by the Institutional Review Board of Naresuan University as all data of our study are deidentified. The analyses were performed at aggregate level (province level) and no confidential Table 1. Variables used to identify risk of infection from GLMM and their descriptions.

Software.

For the maps and for data analysis, we used R program version 4.1.0 [START_REF] Core | R: A language and environment for statistical computing[END_REF] with package tidyverse 1.3.0 45 , stringr 1.4.0 46 , spdep 1.1-5 [START_REF] Bivand | Comparing implementations of global and local indicators of spatial association[END_REF] , sf 0.9-6 [START_REF] Pebesma | Simple features for R: standardized support for spatial vector data[END_REF] , and ggpubr 0.4.0 [START_REF] Alboukadel | ggplot2' Based Publication Ready Plots[END_REF] . The cross-correlation plots were created using R package ggcorrplot 0.1.3 [START_REF] Alboukadel | ggcorrplot: Visualization of a Correlation Matrix using 'ggplot2[END_REF] . For the heat map created using Microsoft Excel version 16.57 (https:// www. micro soft. com/ th-th/ micro soft-365/ excel).

Results

Descriptive results. From 2012 to 2018, there were 20,459 reported leptospirosis cases in Thailand. The annual incidence per 100,000 population (incidence rate) is plotted in Fig. 1. The incidence rate ranged from 3.1 in 2015 to 6.43 in 2012. After 2012, the incidence rate declined till 2014, where it levelled off, then slightly increased after 2016. In addition, the monthly incidence rate was also analysed (Fig. 1). We found that the incidence rates were higher, than the average of 2.47, between June and November, with the highest value in October. Spatial variation was observed in the northern and southern parts Thailand, while the incidence rate in the north-eastern part was observed all year round (Figs. S1 andS2). The incidence rate in the central part was found to be lower than other regions.

Spatial-temporal analysis.

To investigate the spatial autocorrelation, the Global Moran's I values were calculated (Tables S1 andS2). The results showed significant positive values, demonstrating a non-random distribution throughout the period studied, for both annual and monthly incidence rates. To identify high-risk provinces, results of the analysis of local indicators of spatial association were mapped in Figs. 2A and3A. For the annual data, the LISA maps show that high-high clusters appeared in the lower north-east, except in 2012, the year with the lowest Moran's I value. In the south, the high-risk provinces display different patterns. For monthly data, we found high-risk provinces in north-east during rainy season correspond to Fig. 4. The highrisk areas were observed in the south-west correspond to high rainfall (Fig. 4). However, for south-east, the high-risk areas were observed all months except in September and October. Low-risk provinces were mainly detected in central Thailand for both annual and monthly data. Whereas the LISA results of the north showed that there were no significant clusters.

Risk factors.

To select explanatory variables in the model for annual reported data, we calculated the crosscorrelation of the scaled annual variables (Fig. S3). Because precipitation, percentage of flooded area, temperature range, and elevation were highly correlated, they could not be selected in the same model. Elevation, the variable with the lowest univariate model AIC (best model) was kept for model selection. The set of best models (with a ΔAIC < 2 to the model with the lowest AIC) is shown in Table S3 and, following the parsimony principle, we selected the model with the lowest number of variables within this set of best models. Thus, for annual reported data, the final model included annual population density, elevation, and primary rice crop arable areas. These three variables were positively associated with the annual reported cases (Table 2). The set of these explanatory variables is mapped and shown in the supplementary data (Figs. S7-S15). Higher population density was generally found to increase the risk of infection. Elevation was also a predictor for leptospirosis. Specifically in the north-east region, we found that the primary rice crop arable areas were associated with a high probability of leptospirosis (Fig. S13). The estimated annual leptospirosis cases in each year are mapped in Fig. 2B. We found high numbers of predicted cases of leptospirosis in all parts of Thailand except for the Central region, where there were relatively few cases. Figure 5 showed that our model can fit to the reported cases. Note that the weather variables were not included in the final model because of the selection criteria. The same procedure was applied to the monthly data. For monthly reported data, the variables for selecting the final model were temporal correlation, pig density, buffalo density, monthly average soil moisture, monthly temperature range, and monthly NDVI-2 percentage (Table S4). The set of selected variables is also mapped in the supplementary data (Figs. S15-S21). The set of models with a ΔAIC < 2 is shown in Tables S4 andS5. We found that the model using temperature range showed a lower AIC compared to the one using the percentage of flooded area. However, the model using flooded area used less variables, as NDVI-2 was not included. The results from both final models are presented in Table 3 and Table S6. For both models, the temporal autocorrelation parameter was included, indicating the positive correlation of the previous month. The monthly average soil moisture was found to be a positive correlation to the cases, implying the greater the soil moisture the more risk of infection. Temperature range was found to be strongly negative correlation to percentage of flooded area (r ≤ -0.7). Also, the monthly NDVI-2 percentage were negatively with the monthly reported cases. With the alternative final models (Table S6), we found that the percentage of flooded area is positively correlated with monthly reported cases. The results indicated that low temperatures together with a high percentage of flooded area could increase infection risk. The estimated monthly leptospirosis cases are mapped in Fig. 3B. Both models can be used to illustrate the average monthly leptospirosis cases as shown in Fig. 6, Figs. S5,S6.

Discussion

The spatial-temporal analysis of all human leptospirosis cases reported in Thailand from 2012 to 2018 allowed us to highlight the main spatial and temporal patterns and to identify key risk factors associated with leptospirosis infection in Thailand. There are few studies of the spatial-temporal patterns of leptospirosis in the country and none to our knowledge using such a large data set. Although the association of climatic and environmental factors on leptospirosis has been studied in many previous works [START_REF] Chadsuthi | The effects of flooding and weather conditions on leptospirosis transmission in Thailand[END_REF][START_REF] Chadsuthi | A remotely sensed flooding indicator associated with cattle and buffalo leptospirosis cases in Thailand 2011-2013[END_REF] , few studies used remotely sensed data to analyze directly the impact of flooding and took into account primary and secondary rice crop arable land, cultivated areas and yields.

Overall, our analysis showed high incidence rates for specific provinces and seasons in Thailand, suggesting that leptospirosis is still a major public health concern as it included in the prevention research program to control disease by Department of Disease Control, Ministry of Public Health. The spatial autocorrelation analysis highlighted the significance of the annual and monthly spatial clustering of the leptospirosis cases. For the annual data, there are limited differences between the years, as the same provinces (provinces from north-eastern Thailand) have higher risk almost every year. Whereas for the monthly data, the high-risk provinces were different between the months for different parts of the country, especially for the north-eastern and southern parts. It is interesting that in October, which has the highest incidence rate, the high-risk provinces were observed only in north-eastern parts. This may be explained by the process of rice cultivation, which involves many activities such as preparing the land, sowing the wet fields and weeding as discussed below. The hot spots in the north-eastern are almost exclusively rural areas. On the other hand, the hot spots in the south may not be related only to rice cultivation but rather to flooding, as an outbreak occurred 2 weeks after flooding in Nakhon Si Thammarat in January 2017 [START_REF] Chutinantakul | Outbreaks of leptospirosis after a flood in Thung Song District, Nakhon Si Thammarat[END_REF] . It also may have a different seasonality component as compared to the north-eastern part [START_REF]The Climate of Thailand[END_REF] and a different landscape. The low-risk provinces were clustered in the central region of Thailand, which has the highest rice productivity in the country [START_REF] Suebpongsang | Commercialisation of rice farming in northeast thailand in White Gold: The Commercialisation of Rice Farming in the Lower Mekong Basin 39-68[END_REF] .

Our analysis highlights that the areas of different rice cultivated types (in-season and off-season) and rice yield should be considered separately to describe the risk of leptospirosis infection. Consistent with our results, rice farming activities have been identified to be an important risk factor [START_REF] Tangkanakul | Risk factors associated with leptospirosis in northeastern Thailand, 1998[END_REF][START_REF] Sahneh | Investigation of risk factors associated with leptospirosis in the North of Iran (2011-2017)[END_REF][START_REF] Maze | Risk factors for human acute leptospirosis in northern Tanzania[END_REF] . In Thailand, rice farming is the predominant occupation and rice farmers usually cultivate rice in two seasons. Our results demonstrated that the primary rice crop arable area (rice grown using rainwater) was positively associated with the annual leptospirosis incidence. The primary rice crop (in-season) is cultivated during the rainy season (May-October), thus increasing activities in the field, such as fertilizing and ploughing rice in wet fields, leading to a higher risk of exposure to Leptospira. In contrast, the secondary rice crop arable area (off-season) was not considered be a risk factor for leptospirosis as it corresponded to rice growing activities during the non-rainy season. However, rice yield, mainly in the central region, was found to be negatively correlated to primary rice crop arable area. The central region has soil conditions more suitable for increased rice production per area and higher water resources allowing for several rice harvest per year. Our results may suggest that the contamination of Leptospira may depend on the soil characteristics [START_REF] Gamage | Understanding leptospirosis eco-epidemiology by environmental DNA metabarcoding of irrigation water from two agro-ecological regions of Sri Lanka[END_REF] . It may also be due to the different practices for rice cultivation in the central region, which has more mechanized processes than farmers in north-eastern region. Another reason may be the highly effective healthcare system in the central region compared to other parts of the country. The leptospirosis infection could be explained by using land use such as rice growing data and landscape such as elevation. Elevation may be a proxy for increased flooding because lower elevations store rainwater. However, elevation was calculated at province level, the finer spatial resolution should further analysis to be more specific calculation. Due to the limitation of the data (annual data for rice yield), further analysis for the monthly data should be evaluated to calculate the risk factors associated with different seasons.

For the monthly data, we found that using temperature data leads to a minor improvement in the model compared to using the percentage of flooded area. The results showed both parameters can be predictor variables for leptospirosis. Our analyses revealed that temperature has a strongly negative correlation to percentage of flooded area. Usually, in tropical climates, low temperatures correlate to cloudy weather, where there is high rainfall and therefore prone to flood [START_REF] Yong | Variability of diurnal temperature range in east Africa during 1921-2010[END_REF] . The temperature range could be a predictor variable, when the percentage of flooded area is not available. However, using temperature range instead of the percentage of flooding area should be only used in similar climatic areas as Thailand where the model was built and where the correlation between these two factors is likely to hold true. More generally, any extrapolation of this model to other climatic regions should be implemented carefully.

Our finding allows the potential risk of leptospirosis infection to be estimated almost in real-time. Leptospira can survive in freshwater 2 where flooding events could increase the number of pathogenic Leptospira 57 . The average soil moisture, or the soil water capacity [START_REF] Abatzoglou | a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015[END_REF] , was identified as a risk factor for monthly data in this work. Leptospira can survive in soils with a moisture content of ≥ 20% [START_REF] Saito | Comparative analysis of Leptospira strains isolated from environmental soil and water in the Philippines and Japan[END_REF] . Increased soil moisture may increase the survival rate in contaminated soil and water [START_REF] Baquero | Spatiotemporal dynamics and risk factors for human Leptospirosis in Brazil[END_REF] . Note that our results were based on the case reports in the surveillance system. This may not be accurate because some mild symptoms or asymptomatic cases are not going to the hospital, resulting in underreporting of leptospirosis cases. A finer spatial and temporal scale should also be conducted, when well-represented data is available. The analysis could only be performed using the data set by the month and year due to the low number of cases. Regardless of these limitations, our study has provided important knowledge on leptospirosis occurrences by characterizing the hot-spots and key risk factors.

Conclusion

In summary, our study highlighted that when studied over a long period, leptospirosis cases demonstrated a spatial-temporal distribution and showed hot-spots and clustered areas. Regression analysis was used to find the risk factors associated with leptospirosis cases. Working in rice field, which have high soil moisture, could increase risk of infection due to occupational behaviour of working in wet fields facilitating the survival of Leptospira bacteria. The temperature range, which is negatively correlated to the percentage of flooded area was found to be negatively associated with cases. Therefore, flood exposure controls should be used to reduce the risk of leptospirosis infection. These results could be used to improve prevention measures and control actions for public health organizations. The investigation of the associated risk factors can contribute to the prevention and early warning of leptospirosis in Thailand.

Figure 1 .

 1 Figure 1. Annual (A) and monthly (B) leptospirosis incidence rates in Thailand (per 100,000 population).

Figure 2 .

 2 Figure 2. The plots of annual spatial clustering of incidence rate as determined by LISA (A) and annual leptospirosis cases estimated by the final GLMM (B). Maps created using R Program version 4.0.3 (https:// www.r-proje ct. org/).

Figure 3 .

 3 Figure 3. Plots of monthly spatial clustering of incidence rate as determined by LISA (A) and monthly leptospirosis cases estimated by the final GLMM (B). Maps created using R Program version 4.0.3 (https:// www.r-proje ct. org/).

Figure 4 .

 4 Figure 4. Heatmaps of monthly precipitation, monthly percentage of flooded area, monthly temperature range, and average of incidence rate for 3 regions. Red highlights the high values. Blue highlights the low value. The 3 regions were mapped in dark grey. The heat map created using Microsoft Excel version 16.57 (https:// www. micro soft. com/ th-th/ micro soft-365/ excel). Maps created using R Program version 4.0.3 (https:// www.r-proje ct. org/).

Figure 5 .

 5 Figure 5. Plots of comparison between average annual reported leptospirosis cases (A) and average annual fitted leptospirosis cases estimated by the final GLMM (B). Maps created using R Program version 4.0.3 (https:// www.r-proje ct. org/).

Figure 6 .

 6 Figure 6. Plots of comparison between average monthly reported leptospirosis cases (A) and average monthly fitted leptospirosis cases estimated by the final GLMM for temperature range parameter (B). Maps created using R Program version 4.0.3 (https:// www.r-proje ct. org/).

  Temporal correlationThe current number of case reports per population of the province compared with the previous year or month http:// www. boe. moph. go. th/ boedb/ surda ta/ disea se. php? ds= 43

	Variable	Description, (unit)	Reference
	Spatial correlation	Correlation of exponential decay of distance between centroids
		Annual human population per area at each province (population/km 2 )
	Population density	Average human population per area at each province across 7 years (population/ km 2 )	http:// www. world pop. org
		Annual cattle population per area at each province (population/km 2 )
	Cattle density	Average cattle population per area at each province across 7 years (population/ km 2 )	http:// ict. dld. go. th
		Annual pig population per area at each province (population/km 2 )
	Pig density	Average pig population per area at each province across 7 years (population/ km 2 )	http:// ict. dld. go. th
		Annual buffalo population per area at each province (population/km 2 )
	Buffalo density	Average buffalo population per area at each province across 7 years (population/ km 2 )	http:// ict. dld. go. th
	Soil moisture	Annual average of monthly soil moisture at each province Monthly soil moisture at each province, averaged across 7 years (mm)	http:// www. clima tolog ylab. org/ terra clima te. html
	Precipitation	Annual average of monthly precipitation at each province Monthly precipitation at each province, averaged across 7 years (mm)	http:// www. clima tolog ylab. org/ terra clima te. html
	Minimum temperature	Annual average of monthly minimum temperature at each province Monthly minimum temperature at each province, averaged across 7 years (°C)	http:// www. clima tolog ylab. org/ terra clima te. html
		Annual average of monthly maximum temperature at each province
	Maximum temperature	Monthly average maximum temperature at each province, averaged across	http:// www. clima tolog ylab. org/ terra clima te. html
		7 years (°C)
		Difference between maximum annual average temperature and minimum
	Temperature range	annual average temperature at each province Difference between maximum monthly average temperature and minimum	http:// www. clima tolog ylab. org/ terra clima te. html
		monthly average temperature at each province (°C)
	Percentage of flooded area	The percentage of pixels with an MNDWI value greater than or equal to zero (250 m, 8 days) for each month and for each year	MOD09A1
	NDVI-1	The percentage of NDVI pixels that are within 0.2-0.5 (250 m, 8 days) for each month and for each year	MOD09A1
	NDVI-2	The percentage of NDVI pixels that are within 0.5-1.0 (250 m, 8 days) for each month and for each year	MOD09A1
	Elevation	Average elevation of each province (90 m spatial resolution)	https:// srtm. csi. cgiar. org/
	Slope	Average angle of inclination at each province (90 m spatial resolution)	Calculated from elevation
	Primary rice crop arable area	Annual arable area for primary rice crop at each province (km 2 )	http:// www. oae. go. th/
	Primary rice cultivated area	Annual cultivated area for primary rice crop at each province (km 2 )	http:// www. oae. go. th/
	Primary rice yield	Annual primary rice yield of each province (1000 kg/km 2 )	http:// www. oae. go. th/
	Secondary rice crop arable area	Annual arable area for secondary rice crop at each province (km 2 )	http:// www. oae. go. th/
	Second rice cultivated area	Annual cultivated area for secondary rice crop at the province level (km 2 )	http:// www. oae. go. th/
	Secondary rice yield	Annual secondary rice yield of each province (1000 kg/km 2 )

http:// www. oae. go. th/ information was involved. All methods were performed in accordance with the relevant guidelines and regulations. The location and date of onset of illness were also gathered. The maps presented in this paper do not identify patients' addresses.

Table 2 .

 2 Results of the final generalized linear mixed model for annual reported cases. OR odds ratio.

	Variables	OR (95% confidence interval) P-value
	Primary rice crop arable area	1.9576 (1.4398-2.6617)	< 0.0001
	Elevation	1.9259 (1.3125-2.8259)	0.0008
	Annual population density	1.2708 (1.0213-1.5812)	0.0316

Table 3 .

 3 Results of the final generalized linear mixed model for monthly reported cases. OR odds ratio.

	Variables

OR (95% confidence interval) P-value

	Monthly temporal correlation	1.1253 (1.0843-1.1678)	< 0.0001
	Monthly soil moisture	1.3206 (1.2257-1.4230)	< 0.0001
	Monthly temperature range	0.7869 (0.7462-0.8297)	< 0.0001
	Monthly percentage of NDVI-2	0.9533 (0.9186-0.9893)	0.0115

Scientific Reports | (2022) 12:5066 | https://doi.org/10.1038/s41598-022-09079-y

Acknowledgements

We thank Faculty of Science, Naresuan University for their support.

Data availability

The leptospirosis dataset used in the current study are available from Bureau of Epidemiology, Department of Disease Control, Ministry of Public Health, Thailand upon reasonable request. The data supporting the findings can be found in the main paper and in Supplementary Information file.

Author contributions

S.C.: conceptualized, participated in its design, performed the analysis, wrote the first draft and wrote the manuscript. J.C.: conceptualized, participated in its design, analysis and wrote the manuscript. K.C.M.: conceptualized, participated in its design and analysis. S.G.: participated in data extraction, interpretation, and analysis. P.W.: participated in its design and analysis. All authors read and approved the final manuscript.

Competing interests

The authors declare no competing interests.