domain: chim sde

text:
Julius-Alexander Nöpel
email: julius-alexander.noepel@tu-dresden.de

Frédéric Ayela
email: frederic.ayela@legi.cnrs.fr

Experimental evidences of radicals production by hydrodynamic cavitation: a short review

Keywords: Hydrodynamic cavitation, Hydroxyl radicals, Dosimetry, Chemiluminescence, Luminol

The presence of induced radicals and of chemical reactions, in liquids processed in ultrasonic reactors producing acoustic cavitation, has been widely demonstrated and is known as sonochemistry. Otherwise, the number of publications related to chemical reactions caused by hydrodynamic cavitation arising in a liquid flowing out, is more limited. Most of these articles have exposed qualitative results and empirical attempts in order to increase a primary weak efficiency. As the physics of collapsing bubbles in a flow is more complex than the monitoring of acoustic bubbles at a fixed frequency, it was largely assumed, until shortly, that hydrodynamic cavitation has inherent limitations preventing from getting relevant radical yields. However, new efforts have been recently performed to evolve from qualitative to quantitative measurements, in order to be able to know whether hydrodynamic cavitation could become a reliable process for wastewater treatment or not. This short review focuses on recent progresses that have made possible experimental evidences of radical production in hydrodynamic cavitation, and on forthcoming orientations.

Introduction

Cavitation is related to the growth and the collapse of vapor bubbles or vapor clouds in a liquid. The vaporization of the liquid and the condensation of the vapor are thermodynamic changes of phase of the fluid, involving absorbed and emitted latent heat respectively. Mechanical effects induced by the collapses are also likely to occur. Cavitation differs from boiling in the fact that cavitation results from a pressure drop below vapor pressure at a fixed temperature, when boiling results from an increase of temperature above the temperature of saturation, at a fixed pressure.

Acoustic cavitation occurs when oscillating vapor bubbles are induced by acoustic excitation, and monitored at ultrasonic frequencies ranging from 20 kHz to 1 MHz. A global rectified diffusion of matter from the surrounding liquid toward the bubble, during the expansion phases, enhances an increase ISSN (electronic) : 1878-1543 https://comptes-rendus.academie-sciences.fr/chimie/ of the average size of the bubble until a critical radius, that depends on the frequency and on the acoustic pressure, and above which the collapse of the bubble occurs. The consequences of the collapse depend on the frequency, on the intensity of the pressure field, on the presence of dissolved gases, and on the temperature. Imploding cavities often involve the formation of high-pressure shock waves, and extreme temperatures and pressure may be reached. This might not necessarily result in the formation of a supercritical fluid inside the bubble, therefore supercritical oxidation is not a typical outcome in imploding cavities. Low frequency cycles produce large sized bubbles, giving way to violent collapses. As a consequence, some uncondensed water molecules are splitted and produce hydroxyl radicals •OH. High frequency cycles induce smallest bubbles that collapse with a weaker intensity, but they generate more events and they can help for diffusion of radicals toward the liquid. Hydroxyl radicals can recombine and produce H 2 O 2 . In contrast, when adding H 2 O 2 in the liquid, reaction rates increase [1]. Depending on the hydrophilic or hydrophobic nature of a possible pollutant, the optimum working frequency should be a balance between low or high acoustic frequency. Sonoluminescence is another consequence of the collapse, that is caused by the ionization of rare gases present in the core of the bubbles. It should be noted that the best conditions to observe sonoluminescence are not similar to those enhancing sonochemical activity.

Hydrodynamic cavitation is the consequence of the Bernoulli law, when a flowing liquid passes through a local constriction. The increase of the dynamic pressure 1/2 • ρu 2 , where ρ is the density of the fluid and u its average velocity, goes with a decrease of the static pressure p, as illustrated in Figure 1. Downstream the constriction, where the velocity of the flow recovers a lower value, the static pressure increases but reaches a smaller level than upstream, because of the singular pressure loss due to the constriction. Inception of cavitation becomes possible if the pressure value falls below a threshold value, which is generally the vapor pressure p vap . Vapor bubbles and vapor clouds, which arise in the low pressure area, collapse downstream where the static pressure has raised above p vap .

Different flow regimes can induce cavitation: around an emerging liquid jet, around a hydrofoil, through an orifice or a diaphragm, through the throat of a venturi. Shear cavitation is induced in low pressure vortices in shear flows; sheet cavitation appears as a vapor pocket emitting bursts of vapor clouds at a shedding frequency that has no relation with the frequency monitored in ultrasonic cavitation (Figure 2).

Hydrodynamic cavitation is simple to execute. However, the physics of the growth and collapse of multiple cavitating bubbles in a hydrodynamic flow is drastically more complex than in ultrasonic reactors. A cavitation number, noted σ, is a dimensionless parameter that can be used to quantify the possibility for cavitation to occur in a fluid flow system. With the notations of Figure 1, σ is defined here by:

σ = p downstream -p sat 0.5ϱu 2
The condition σ < 1 is required for the cavitation to become possible. In the literature, some alternative definitions of the cavitation number may be found, but all agree with the fact that σ can be expressed as a function of the ratio of the output pressure p out divided by the flow rate. The cavitation number only means that a cavitating flow regime is reachable by increasing the flow rate and/or by lowering the output pressure and that no cavitation should occur when σ > 1. In microfluidic systems, we have calculated that the critical cavitation number below which the pressure falls under p sat depends on the size of the constriction [2]. More recently, Sarc et al. have highlighted, that describing the hydrodynamic cavitation by the solely value of σ was inadequate [3]. Moreover, some factors affecting chemical reactions, such as turbulence, dissolution of gas, mass transfer rates, are not relevant with the cavitation number. So, it is suggested [3] to consider several parameters, in order to address the issue of different experiments.

We consider that the cavitation number, the maximum average velocity and the hydraulic power, that is the product of the total pressure drop (in J/m 3 ) by the volumetric flow rate of the liquid phase (in m 3 /s), are the quantities that allow comparison between hydrodynamic flows. But the frequency, which is a crucial parameter in acoustic cavitation, has no obvious peer in hydrodynamic cavitation. Such a lack of a parameter able to evaluate an average shape or size of the bubbles may result in the fact that the experimental attempts to study chemical effects induced by hydrodynamic cavitation could not be based on a firm model.

Chemistry induced by hydrodynamic cavitation

Twenty-five years ago, Suslick and co-workers underlined that reports on chemistry induced by hydrodynamic cavitation had been extremely limited [4]. They tested the Weissler reaction, as a dosimeter for chemical reactions, with a solution passed through a high pressure microfluidizer in which cavitation was likely to occur, with working pressures running from 100 bars to 1300 bars. The tri iodide formation rate increased linearly with the upstream pressure. The respective influence of dissolved gases and liquid temperature were also considered, and led to evolutions similar as those observed in acoustic cavitation. That was the first firm demonstration of chemical similarities between acoustic and hydrodynamic cavitation. Unfortunately, the flow rates were not expressed in the paper, and it is not possible to evaluate the hydraulic powers involved in these experiments.

Later, Kalumuck and Chahine tested the degradation of p-nitrophenol aqueous solutions in both acoustic and hydrodynamic reactors [5]. They compared the power of the hydraulic set up to the calorimetric power of the acoustic reactor, and concluded to a better efficiency of the hydraulic reactor. Temperature and pH effects highlighted the same trends, in acoustic and hydrodynamic cavitating regimes. The cavitating jet was monitored up to 65 bars. Comparing the application of cavitation reactors for water disinfection, Gogate [6] enhanced an empirical arrangement of multiholes on a plate, taking into account the relative free area offered to the flow. General conclusions were not supported by quantitative values of parameters such as inlet pressure or average maximum velocity. Perhaps in order to treat flow rates higher than what is reachable with high pressure fluidizers, the subsequent efforts, which were consacred, used mild experimental conditions. Wang et al. studied the degradation of rhodamine B using swirling jet-induced cavitation, at an inlet pressure of 6 bars and a hydraulic power of around 600 W [7]. Brauetigam et al. used a design inspired by the multiholes scheme [8,9] working on a very limited inlet pressure range (from 1.4 bars to 2.6 bars [8], and then up to 10 bars [9]), with hydraulic powers between 60 W and 80 W. They noticed that the centrifugal pump present in the loop could be a source of cavitation, in competition with the reactor. Capocelli et al. undertook simulation and experimental comparison of chemical effects of hydrodynamic cavitation in a Venturi reactor [10]. They conducted degradation tests on p-nitrophenol between 2 bars and 6 bars, and observed a maximum efficiency around 4.5 bars. That optimum was correlated to similar trends obtained by other groups, but with other reactors and other pollutants. Zupanc et al. performed a broad study of shear-induced hydrodynamic cavitation as a tool for pharmaceutical pollutants removal [11]. They stated that the amount of free radicals is related to the aggressiveness of the cavitation, and that shear induced cavitation is more aggressive than sheet cavitation that is formed in the throttle of a Venturi. That conclusion is not shared by Rajoriya et al. [12] who have observed a higher degradation rate of orange-G dye when using slit Venturi at an optimum pressure of 3 bar. The authors state that cavitation in a Venturi profile results into a large magnitude collapse, as long as a choked cavitation regime is not reached. The authors present also an overview of experimental results devoted to waste water treatment by hydrodynamic cavitation through Venturi and orifice nozzles.

A larger spectrum of references may also be found in a short review article published by Tao et al. in the same year of 2016 [13], and more recently by Das et al. [14]. The authors review the possible reactors geometries and list the main operating parameters that are inlet pressure, bulk liquid temperature, effect of pH and of initial pollutant concentration, effect of the initial radius of the nuclei, effect of the surface tension. They agree that qualitative results have been obtained, but they underline the lack of quantitative analysis. Each research team has developed a selfmade reactor used to degrade a given pollutant, and it is difficult to draw general conclusions for such a set of data. Similar conclusions have been drawn in most recent review articles [15,16]. Ciriminna et al. stress that the chemical effects of hydrodynamic cavitation are the consequence of hydroxyl radical formation, and that a direct measurement of such a concentration is hardly possible due to their high reactivity [15]. Hydroxyl radicals are short-lived oxidants that degrade within nanoseconds [START_REF] Dorfmann | Reactivity of the hydroxyl radical in aqueous solutions[END_REF]. Therefore, it is difficult to detect them directly. The dosimetry of hydroxyl radicals comes from indirect measurement, funded on chemical reactions involving new products, or on the emission of light.

Dosimetry of hydroxyl radicals produced by hydrodynamic cavitation

Different dosimetry methods have been developed in sonochemistry in order to quantify the chemical activity and the amount of hydroxyls radicals produced by cavitation. Most of them have then been tested in hydrodynamic cavitating flows. The reliability of these techniques, that can only detect radicals which have not been involved in recombination, is still an open question in sonochemistry [START_REF] Rajamma | [END_REF]. The data concerning dosimetry of hydrodynamic cavitating flows are equally scattered. Iodide dosimetry has been largely used as a test reaction. It has been inspired by the Weissler reaction, where the yield of iodine is proportional to the time of ultrasonic irradiation. When potassium iodide KI is initially present in the solution submitted to cavitation, iodide can scavenge the hydroxyl radicals •OH and form the tri iodide I 3 -, which is detected by UV visible absorbance at λ = 353 nm [START_REF] Rooze | Ultrason[END_REF]20]. Otherwise, a part of the processed solution may be withdrawn and mixed with a KI solution. The oxidation of iodide to tri iodide I 3 -is then caused by hydroxide peroxide H 2 O 2 resulting from the recombination of hydroxyl radicals [21]. But iodide dosimetry is non-specific to •OH radicals, since iodide ions may be oxidised by other radicals and by hydrogen peroxide. The presence of saturated oxygen in the cavitating liquid is likely to help the formation of peroxide hydrogen, and to overestimate the amount of hydroxyl radicals generated by cavitation. Moreover, in a global study devoted to compare the respective efficiencies of acoustic and hydrodynamic cavitation, Morison and Hutchinson enlightened the presence of a consumption reaction in hydrodynamic flows, leading to a limit of the I 3 concentration [20]. That could explain why experiments of Suslick et al. [4] could detect tri-iodide only above very high input pressures. Such a consumption of tri iodide ions occurred even if the flow was driven below the onset of cavitation, but was not present in acoustic cavitation. Their conclusion stated that the Weissler reaction is not reliable in order to know the effectiveness of hydrodynamic cavitation as an advanced oxidizing process.

The Fricke dosimetry is funded on the oxidation of Fe 2+ to Fe 3+ [22]. As for the Weissler reaction, other species than •OH radicals may be responsible for the oxidation reaction. Terephtalic acid TA has been used as an •OH radical scavenger [23,24]. The reaction forms hydroxyterephtalic acid (HTA), which is strongly fluorescent and stable for hours. Terephtalic acid dosimetry is presented as specific for hydroxyl radicals [23]. Matthews, working on the radiation chemistry of the terephtalic dosimeter, observed that the yield of HTA is affected by the oxygen concentration [25]. Oxygenation increases HTA, but in a solution sparged with nitrogen gas, the yield of HTA falls. In acoustic cavitation, the fact that a very low •OH yield is obtained by this dosimetry method is still an open question [START_REF] Rajamma | [END_REF]. In hydrodynamic cavitation, Wang et al. [7] used TA dosimetry to conclude to a synergetic effect between cavitation and H 2 O 2 for the degradation of rhodamine B, with fluid pressures ranging from 2 to 6 bars.

Arrojo et al. proposed the application of salicylic acid dosimetry to characterize both ultrasonic and hydrodynamic cavitation [23]. The reaction with •OH produces mainly 2,5 DHB (dihydroxybenzoic acid), that is fluorescent. The reactor was a venturi, developing a hydraulic power around 500 W at a maximum inlet pressure of 4.7 bars. The authors underlined that as salicylic acid is hydrophobic, molecules should migrate toward the bubbles interface. That would enhance the trapping efficiency, as the life time of •OH is very short. A new question then arises, that is to know the ratio between the real number of hydroxyl radicals initially present inside collapsed bubbles, and the effective number of radicals able to react with a hydrophobic or with a hydrophilic molecule present in the liquid phase. Salicylic acid dosimetry was also used by Amin et al. [26], who tested different shapes of orifices (circular, ellipsoidal, rectangular). But unfortunately, there is a lack of hydrodynamic data (such as pressure drop/flow rate relationships, critical pressure drop above which the onset of cavitation occurs), that cannot be related to chemical activity and to the related measurements. More recently, Zupanc and co-workers [27] used salicylic acid dosimetry with a Venturi equipped with a rectangular slit. They highlighted the fact that an excessive concentration of salicylic acid should negatively affect the surface tension and, as a matter of fact, the violence of the collapse. Care must be taken, so that the scavenger should not affect the radical production.

Coumarin dosimetry has been recently used in cavitation. After an exposition to ultrasound waves at 500 kHz, an aqueous coumarin solution has produced hydroxylated products and has presented a decrease of the absorbance of coumarin, that are the consequence of reactions with hydroxyl radicals [28]. Among these products, umbelliferone (also named hydroxylated 7-hydroxycoumarin, 7-OHC) is fluorescent. Coumarin was then used by Srinivas et al. [29] and by De-Nasri et al. [30] to quantify the •OH radical generation in hydrodynamic cavitating vortex diode reactors. Coumarin was also used for the detection of radicals in irradiated solutions [31]. Once again, it is not obvious to know the exact relationship between the amount of •OH molecules involved in the reaction producing the fluorescent 7-OHC component, and the total number of •OH radicals.

A variety of other methods have been published. For example, the transformation of benzene into phenol due to the reaction with hydroxyl radicals can be used as a dosimeter [32]. The correlation between dye degradation experiments, with Congo red used as a reference substance, and chemiluminescence of luminol was used to evaluate degradation effects and cavitation intensity [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF].

Chemiluminescence of luminol

The above dosimetry methods require a regular extraction of the liquid under test, and do not provide any information neither about the fluid mechanism conditions nor about the localization of the chemical activity. An indirect determination of the presence of •OH hydroxyl radicals and an in-situ visualisation of their localization, may be performed with luminol aqueous solutions. It is known for a long time that aqueous solutions of luminol emit light when exposed to acoustic waves, which are strong enough to generate cavitation, and that such a chemiluminescence results from reactions with hydroxyl radicals [START_REF] Taylor | [END_REF]. Chemiluminescence of luminol has also been enhanced in single bubble sonoluminescence experiments [35]. The method of chemiluminescence of luminol gives a statement about the spatial distribution of the reaction area and is suitable as a marker for the spatial characterization of the hydroxyl radical concentration, which is also quantifiable for measurements. The observation of the intensity of the radical production becomes possible inside devices equipped with optically transparent walls, because when cavitation occurs and hydroxyl radicals are formed, the decomposition of luminol emits visible blue light above the UV-spectrum.

Surprisingly, using luminol as the working fluid in hydrodynamic cavitating experiments has been firstly presented only six years ago [36]. Gathering relatively high volumes of solution together with transparent reactors able to withstand high pressure levels, is perhaps the reason why there has been such a delay between acoustic and hydrodynamic cavitation experiments using luminol. The measurement of chemiluminescence of luminol is now a well-established proof for hydroxyl radical production by hydrodynamic cavitation even if the number of papers dealing with luminol cavitating flows remains up to now limited. On the one hand, several experiments have been performed under moderate pressure drops through microchannels (hydrodynamic cavitation "on a chip") [37][38][39] and mini channels [40]; on the other hand, another set of experiments at macroscale and under inlet pressures up to 40 bars have been published [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF][START_REF] Nöpel | Proceedings of the 5th International Conference on Experimental Fluid Mechanics (ICEFM)[END_REF][START_REF] Deggelmann | [END_REF].

The existence of luminescent scavengers of hydroxyl radicals makes possible to count emitted photons directly. However, the relationship between the number of radicals reacting with luminol and the number of collected photons is not obvious. Chemistry of the oxidation of luminol is complex and the emission of light requires both the presence of hydroxyl radicals together with superoxide anion O 2 -• [43,44]. As a matter of fact, the deficit of superoxide O 2

•-, even if •OH are still present, will stop the luminescence, and one should underestimate the number of •OH radicals. Moreover, when using a photomultiplier tube (PMT) as a photon detector, one has to take into account the correction due to the collection solid angle, the global quantum efficiency of the photocathode of the PMT for the chemiluminescent spectrum, the quantum yield of the luminol chemiluminescence reaction. Whatever the size of the reactor is, the test rig must be completely darkened in order to capture the relatively weakly emitted light. The knowledge of a number of radicals detected during one second, and the knowledge of the hydraulic power required to drive the flow make possible the calculation of the so-called G-value, that is the number of •OH moles produced per Joule of consumed energy. By working with hybrid silicon-Pyrex micro channels with a micro diaphragm inside, Podbevsek et al. calculated a minimum radical production of ≈ 6 × 10 6 •OH/s for a hydraulic power of ≈ 0.4 W, that corresponds to G ≈ 2.50 × 10 -17 mole/J [37].

Similar G-values were also recorded with another opto-electronic set up and with other micro reactors (microdiaphragms and micro steps profiles) developing shear rate cavitation [38]. Recently, the chemiluminescence emitted from a larger set of microand mini hydrocavitating reactors has been investigated [40]. The devices presented semi convergentdivergent constrictions, with hydraulic diameters ranging from 100 µm to 540 µm. For the largest diameter, the maximum recorded flow rate was 50 L/h under a driving pressure of 10 bars. The flow dynamics was rather different from that inside micromachined silicon-Pyrex devices, because cloud shed-ding was present at frequencies ranging from 150 Hz to 1400 Hz. It is noteworthy to find that the G-values calculated from the reactive oxidising species production were similar (3.6 × 10 -17 moles/J < G < 6 × 10 -17 moles/J) to those calculated from former experiments. Moreover, the radical production exhibits a linear evolution as a function of the hydraulic power, and no scale effect was noticeable: the G values of the microreactors were roughly similar to the G values of the minireactors. For a small as for a large device, the G-value was the highest for the smallest constriction under test. To our knowledge, only Arrojo et al. had up to now calculated a G-value from hydrocavitation experiments, and found G ≈ 10 -12 mole/J with salicylic acid as a scavenger used in a Venturi [23]. There are five orders of magnitude between these respective G-values. The microfluidics experiments of Podbevsek et al. [37,39] and of Perrin et al. [38] were mainly devoted to demonstrate the effective presence of hydroxyl radicals inside cavitating flows, and they were performed just above the critical flow rate corresponding to the inception of cavitation, when the cavitating flow was not fully developed. But the recent results obtained at larger scale confirm the former G values [40]. Possible explanations are as follows: the value of the quantum yield of luminol is under evaluated; the luminescence is limited by the superoxide anion production; the hydrophobic salicylic scavenger react with a larger number of radicals, which have disappeared when using luminol. Anyway, complementary experiments at micro and miniscale, with complementary geometries, are scheduled in order to draw firm findings.

Even if chemiluminescence of luminol cannot provide the real number of hydroxyl radicals, it makes possible a mapping of the chemical activity together with a mapping of the two-phase liquid-vapor flow. Supporting the fact that cavitation occurs inside a "lab on a chip", Podbevsek et al. used a confocal microscope set up to draw a local mapping of the chemiluminescent intensity, enhancing local gradients of intensity. Such mappings could be compared to void fraction mappings, also obtained from confocal microscopy with fluorescent hydrophilic nanotracers [39]. An example of a chemiluminescent mapping is presented in Figure 3, with a snapshot of the corresponding two-phase flow. The main conclusion was that the area where the maximum chemical activity was recorded (due to the presence A correlation between the volume of chemiluminescence due to luminol and the volume of cavitation has been studied by Deggelmann et al. [START_REF] Deggelmann | [END_REF] and Nöpel et al. [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF][START_REF] Nöpel | Proceedings of the 5th International Conference on Experimental Fluid Mechanics (ICEFM)[END_REF], but on a larger scale. The authors applied the method of chemiluminescence of luminol, which occurs in the visible range of light due to the reaction with hydroxyl radicals. The emitted light was then recorded with a digital SLR camera.

The investigations about the luminescent area and the cavitating vapor cloud were performed initially in a cylindrical reactor with an inner diameter of 15 mm and a length of 100 mm [START_REF] Nöpel | Proceedings of the 5th International Conference on Experimental Fluid Mechanics (ICEFM)[END_REF][START_REF] Deggelmann | [END_REF]. The flow was accelerated through a circular orifice with diameter sizes of d = 1 mm or d = 1.7 mm. In [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF] same diameters were adapted, but the figuration changed to a rectangular reactor with a cross-sectional area of 15 × 15 mm 2 , leading to less light distortion due to curved surfaces. The measurement configuration of recording the chemiluminescence consisted of a digital SLR camera and a light sensitive lens, which was aligned orthogonally to the reactor walls.

The experimental setup with the reactor was shielded from external light influence. In this darkened room, blue light emission from the chemiluminescence of luminol was recorded by long exposure times of 5 min [START_REF] Nöpel | Proceedings of the 5th International Conference on Experimental Fluid Mechanics (ICEFM)[END_REF] to 10 min [START_REF] Deggelmann | [END_REF]. This method is relatively inexpensive and straightforward to perform. The solution here consisted of 2 g/l luminol and 7.5 g/l sodium hydrogen carbonate, which are both dissolved in deionised water by stirring. The solution has a pH value of pH = 9.3 [START_REF] Nöpel | Proceedings of the 5th International Conference on Experimental Fluid Mechanics (ICEFM)[END_REF].

The method of measuring the chemiluminescence of luminol involved the digital single lens camera recording an integral value of the blue light emission over the exposure time and space. Figure 4 shows the distribution of vapor bubbles (top) and the distribution of light emission of chemiluminescence by luminol reacting with hydroxyl radicals [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF]. In Nöpel et al. [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF][START_REF] Nöpel | Proceedings of the 5th International Conference on Experimental Fluid Mechanics (ICEFM)[END_REF] the exposure time was set to 5 min without additional substances. A camera with a light sensitive lens (Canon EOS 650D, EF 50 mm 1:1.8 STM) was used for this purpose. A different camera and lens (Canon EOS 1100D, EF-S 18-55 mm IS II) was used in [START_REF] Deggelmann | [END_REF] and exposure time was set to 10 min resulting in different area sizes. Also a low amount V = 1 ml of hydrogen peroxide solution (30%) of was added in [START_REF] Deggelmann | [END_REF], which enhanced the reaction. The changes towards smaller areas were noted in [START_REF] Deggelmann | [END_REF]. This was due to the different measurement equipment.

The original images are first brightened and binarized, all similarly, using a threshold for evaluation. From the binary images, the pixel area of Original image is brightened [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF].

the blue light emission can be determined, which represents two-dimensional information about the reaction area. Based on this information, the length extent of the area is used to calculate an approximate volume in which reactions are assumed to be occurring [START_REF] Nöpel | Proceedings of the 5th International Conference on Experimental Fluid Mechanics (ICEFM)[END_REF].

Figure 5 displays the dependency of the reaction volume marked by chemiluminescence over hydraulic power input. The inception of cavitation occurred at ∆p = 7 bar for d = 1 mm. For pressure differences below 10 bar, no light was recorded since the cavitation area was short and vapour density optically low. From 10 bar to 30 bar the optically visible volume of chemiluminescence increases nearly linear.

The higher the pressure difference ∆p or the hydraulic power P, the greater the cavitation volume of gas bubbles. A larger volume of cavitating bubbles leads to a higher implosion rate, which must tend to produce more radicals. However, the pressure difference increases the fraction of vapour and air bubbles. This may result in differences in the spatial extent of the recorded area of chemiluminescence due to scattering and reflexion of light by the gas fraction. This causes a certain limitation for the application of measuring the strength of cavitation for degradation. In [START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF] the method is used to proof the occurrence of oxidation as the leading mechanism of degradation.

Conclusion

The production of hydroxyl radicals in cavitating flows has been attested by an amount of experiments, each of them performed with a specific hydraulic set up and with a fixed chemical scavenger. Compared with ultrasonic cavitation, that involves a relatively limited number of parameters, hydrodynamic cavitation may be monitored in a lot of kinds of devices, under different complex flow regimes. The number of publications related to chemical activity in hydrodynamic flows has increased these past last years, because possible applications to wastewater treatments are timely. Concerning the hydroxyl radical production, the pros and cons of related scavengers have been now assessed. Recent progresses with luminescent scavengers have made possible live and in-situ observations, facilitating coupled studies between chemistry and fluid mechanics. If the rare quantitative data published up to now agree with the fact that hydrodynamic cavitation could display a low radical efficiency, it is nevertheless believed that any further research programs should develop an experimental, structured and multidisciplinary strategy, in order to reach a full understanding of the physics inside collective collapsing bubbles. For each kind of hydraulic device, the largest number of scavengers should be used and compared, for different cavitating flow regimes. Reciprocally, each scavenger should be tested in different devices, at different scales, and submitted to shear and to sheet cavitation. Ideally, every experiment could be managed in order to rise a quantifiable radical efficiency, expressed in •OH moles/J, that would be a common reference for any comparison.

Figure 1 .

 1 Figure 1. Sketch of the pressure line evolution along a liquid flow encountering a constriction.

Figure 2 .

 2 Figure 2. Snapshots of shear (up) and sheet (down) cavitating flows.

Figure 3 .

 3 Figure 3. Local repartition of emitted photons per second (up) compared to a snapshot of the corresponding cavitating flow (down) recorded on a microdiaphragm. The fact that the coloured area oversizes the channel area is the consequence of the diffusion of light, and of the lack of filters removed in order to increase the sensitivity.

Figure 4 .

 4 Figure 4. Inside the cavitation reactor with nozzle d = 1 mm at ∆p = 30 bar. Flow from left to right. Top: Time averaged spatial distribution of vapor bubbles representing area of cavitation. Acquired by side illumination with high power lighting and a camera with exposure time of 11 µs. Bottom: Time averaged spatial distribution of chemiluminescence of luminol (blue light) representing the reaction area.Original image is brightened[START_REF] Nöpel | 11th International Symposium on Cavitation[END_REF].

Figure 5 .

 5 Figure 5. Volume of chemiluminescent part emerging from orifice d = 1 mm normalized with reactor volume over hydraulic power P in W starting from ∆p = 0 bar up to 30 bar in steps of 10 bar [41].

Acknowledgments

The authors gratefully acknowledge the support provided by their respective laboratories. Considerable appreciation shall be expressed to Professor J. Fröhlich, Director of the institute and Head of the chair of fluid mechanics and Dr. F. Rüdiger, Head of the group of experimental fluid mechanics, who provided the opportunity for experimental research on cavitating multiphase flows at the institute of fluid mechanics of the TU Dresden. FA, from the Laboratoire des Ecoulements Géophysiques & Industriels, expresses his sincerest appreciation to Dr. L. Perrin, Dr. D. Colombet, Dr. D. Podbevsek and Dr. G. Ledoux for the useful help and effective collaborations that have led to the success of sophisticated experiments devoted to hydrodynamic cavitation on a chip.

Conflicts of interest

Authors have no conflict of interest to declare.


Fabrizio Falasca
email: fabrifalasca@gatech.edu

Julien Crétat

Pascale Braconnot

Annalisa Bracco

Spatiotemporal complexity and time-dependent networks in mid-to late Holocene simulations

In climate science regime transitions include abrupt changes in modes of climate variability and shifts in the connectivity of the whole system. While important, their identification remains challenging. This paper proposes a new framework to investigate regime transitions and connectivity patterns in spatiotemporal climate fields. Firstly, local regime shifts are quantified by means of information entropy; secondly, their spatial heterogeneity is examined by identifying the underlying spatial domains of the entropy field; finally, a weighted, direct and time-dependent network is inferred to capture the linkages between these domains. The spatiotemporal variability in sea surface temperature (SST) in two simulations of the last 6000 years is investigated with the proposed approach. The largest regional regime shifts emerge as abrupt transitions from low to high-frequency SST oscillations, or vice versa, in both simulations. Furthermore, the variability in time of the two climate networks is studied in terms of their network density. Generally, rapid and sudden transitions in the degree of connectivity of the system are observed in both simulations but, in most cases, at different times, with few exceptions. This suggests that our ability to predict the climate system may be hampered by its inherent complexity resulting from internal variability.

Introduction

Our ability to predict Earth's climate is limited by the nonlinear and intrinsically complex interactions between its components. Quantifiable changes in the system dynamics, i.e. dynamical transitions or climate regime shifts, are manifestations of such complexity. A global-scale example of major climate shifts are the Snowball Earth glaciations that may have occurred 650-700 million years ago [START_REF] Hoffman | The snowball earth hypothesis: testing the limits of global change[END_REF][START_REF] Pierrehumbert | Climate of the neoproterozoic[END_REF][START_REF] Lucarini | Edge states in the climate system: exploring global instabilities and critical transitions[END_REF][START_REF] Lucarini | Transitions across melancholia states in a climate model: reconciling the deterministic and stochastic points of view[END_REF]. In addition, the climate system has experienced a large variety of regional regime shifts, such as the Dansgaard-Oeschger events [START_REF] Dansgaard | A new greenland deep ice core[END_REF][START_REF] Johnsen | Irregular glacial interstadials recorded in a new greenland ice core[END_REF] or the abrupt hydrological fluctuations in the African and Asian monsoon regions [START_REF] Gasse | Hydrological changes in the african tropics since the last glacial maximum[END_REF][START_REF] Morrill | A synthesis of abrupt changes in the asian summer monsoon since the last deglaciation[END_REF], the sudden cooling of the North Pacific sea surface temperature in the mid-Holocene [START_REF] Mix | Rapid climate oscillations in the northeast pacific during the last deglaciation reflect northern and southern hemisphere sources[END_REF] or the 4.2-kiloyear BP aridification event during the Holocene epoch [START_REF]Cultural responses to climate change during the late holocene[END_REF] to cite a few 1 . In recent years, methods stemming from non-equilibrium statistical mechanics, information theory, and network analysis have proven valuable for understanding and investigate complex climate transitions [START_REF] Balasis | Statistical mechanics and information-theoretic perspectives on complexity in the earth system[END_REF][START_REF] Donner | Recurrence-based time series analysis by means of complex network methods[END_REF]. Examples include the work by Donges et al. [START_REF] Donges | Nonlinear detection of paleoclimate-variability transitions possibly related to human evolution[END_REF], who investigated climate fluctuations in Northern and Eastern Africa and linked them to rapid changes in hominid evolution, or by Donges et al. [START_REF] Donges | Non-linear regime shifts in holocene asian monsoon variability: potential impacts on cultural change and migratory patterns[END_REF], who showed that nonlinear regime shifts of the Asian monsoon were associated to the migration and collapse of ancient societies. In the study of dynamical systems, the permutation entropy metric [START_REF] Bandt | Permutation entropy: a natural complexity measure for time serie[END_REF] is often used in nonlinear time series analysis to quantify the complexity of a dynamical system, with changes in complexity being linked to regime shifts. This metric has been applied recently to climate fields to studying the complexity of a proxy record of the El Niño Southern Oscillation (ENSO) [START_REF] Saco | Entropy analysis of the dynamics of el niño/southern oscillation during the holocene[END_REF] and to identifying anomalies and post-processing issues of deep polar ice-cores [START_REF] Garland | Anomaly detection in paleoclimate records using permutation entropy[END_REF]. An alternative, more traditional measure of complexity in dynamical systems is provided by Lyapunov exponents [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF] that quantify the rate of separation of infinitesimally close trajectories. However, recently developed entropy quantifiers are easier to estimate and preferable when dealing with arbitrary real-world observations and large datasets [START_REF] Bandt | Permutation entropy: a natural complexity measure for time serie[END_REF]. Independently of the methodology used, the majority of the studies employing nonlinear data analysis techniques for studying climate dynamics have focused so far on single time series (i.e., on paleoclimate proxies or their modeled equivalent). Here, we propose a novel approach that allows for extending nonlinear time series analysis to the investigation of dynamical transitions in spatiotemporal fields. The methodology couples an entropy quantifier recently proposed by Corso et al. [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF] with δ-MAPS a dimensionality reduction scheme developed by Fountalis et al. [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF]. Given a spatiotemporal field, our framework allows to quantify the local (i.e., for each grid cell) time evolution in information entropy and identify spatially contiguous sets of grid cells (domains) with highly correlated entropy variability. Changes in the information entropy of a time series reflect changes in its recurrent structure and hence in its system's dynamics. It follows that the identified domains consist of grid cells sharing the same dynamical transitions. Studying the average entropy signal inside domains enables to assess the timing, strength and location of climate regime shifts. In this paper we generally focus on abrupt regime shifts, broadly defined as fast and sudden transitions in the system's dynamics. A strong abrupt regime shift in the sea surface temperature (SST) field, and reported in this work, would consist in a sudden transition in regional dynamics from high-to low-frequency variability, or vice versa, that occurs over a period much shorter than the length of the integrations, for example over 200 to 300 years. Additionally, the dynamics of the climate system establish linkages between domains, and complex network science provides a well-established methodology to quantify these connections [START_REF] Donner | Understanding the earth as a complex system -recent advances in data analysis and modelling in earth sciences[END_REF][START_REF] Donges | Complex networks in climate dynamics[END_REF]. Here we adopt a functional network inference that allows to quantify the magnitude of connections between domains, their directionality, and their evolution in time [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF]. In our framework, domains are considered as static nodes of the climate network and links quantify their connections [START_REF] Barabási | Network science[END_REF][START_REF] Newman | Networks: An introduction[END_REF]. The strength of our approach relies in exploring climate variability simultaneously in space and time.

To showcase this novel methodology, we analyze two global ocean-atmosphere coupled simulations of the last 6000 years performed with different versions of the IPSL Earth system model [START_REF] Dufresne | Climate change projections using the IPSL-CM5 earth system model: From CMIP3 to CMIP5[END_REF][START_REF] Braconnot | Strengths and challenges for transient mid-to late holocene simulations with dynamical vegetation[END_REF]. These simulations test the response of the climate system to the slow variations in Earth's orbit and to changes in trace gases. They produce a wide range of variability resulting from the combination of the response to external forcing (Earth's orbit and gases) and intrinsic climate variability [START_REF] Braconnot | Strengths and challenges for transient mid-to late holocene simulations with dynamical vegetation[END_REF][START_REF] Braconnot | Impact of multiscale variability on last 6000 years indian and west african monsoon rain[END_REF]. These two long, global climate simulations provide an ideal testbed for assessing the efficiency of the new method and its ability in identifying major shifts in the simulations.

The paper is organized as follows: Section 2 introduces the proposed approach; the climate simulations analyzed are described in Section 3; results are presented in Section 4 and a comparative application to a reanalysis data set follows in Section 5. A discussion of the method and its preliminary applications concludes the work.

Methodology and strategy

The proposed strategy builds upon concepts of recurrence plots, information entropy and δ-MAPS.

Here we briefly introduce these tools and present the rationale behind and strategy followed by our comprehensive approach. More details on each component can be found in Marwan et al. [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF], Corso et al. [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF] and Fountalis et al. [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF].

Recurrence plots and information entropy

A recurrence plot (RP) is a nonlinear time series analysis method introduced by Eckmann et al. [START_REF] Eckmann | Recurrence plots for the analysis of complex systems[END_REF] that allows for visualizing a d -dimensional dynamical system in terms of its recurrences, a fundamental property of these systems [START_REF] Poincaré | Sur le probleme des trois corps et les équations de la dynamique[END_REF][START_REF] Cvitanović | Chaos: Classical and quantum[END_REF]. Given a trajectory x i of a dynamical system in a d -dimensional state space at time i, a RP is an N × N matrix of 1 and 0 given by:

RP i,j ( ) = Θ( -x i -x j ), x i ∈ R d , i, j ∈ [1, N ] . (1)
Here is the threshold distance and defines the neighborhood of a state x i , Θ is the Heaviside function, • is a norm and N is the number of states considered. The analysis of structures (such as diagonal, vertical or horizontal lines) in a RP is known as recurrence quantification analysis (RQA) and has found numerous applications [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF]. RQA allows for the definition of different measures of complexity of time series. Here we focus on one such measure, developed by Corso et al. [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF] and based on the probability of occurrence of microstates in a RP. A microstate of size M is defined as an M × M matrix inside the RP. The total number of configurations of 1 and 0 in a microstate of size M is M * = 2 M 2 . It is then possible to define a probability of occurrence P k of a microstate k as P k = n k M * , with n k being the number of occurrences of the k-th microstate in the RP. The information entropy of the time series is then given by:

S(M * ) = - M * k=1 P k log P k . (2)
The maximum entropy corresponds to the case in which all microstates appear with the same probability P k = 1 M * . It can be analytically calculated as S max (M * ) = M log 2 and used to normalize the entropy, so that S = 0(1) implies total predictability (unpredictability) of the system dynamics. The number of admissible microstates M * grows exponentially as a function of M , but Corso et al. [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF], showed that in practice only few of them actually populate the RP. It is therefore reasonable to sample M of the possible microstates in the RP and to expect a rapid convergence to S(M * ). In this case, the maximum entropy will be equal to S max ( M ) = log M . The computed information entropy is a measure of the complexity of the time series examined and from now on the two terms will be used interchangeably. There are other methods to estimate complexity of time series, for example related to RQA [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF], recurrence network analysis [START_REF] Donner | The geometry of chaotic dynamics -a complex network perspective[END_REF] and information theory [START_REF] Balasis | Statistical mechanics and information-theoretic perspectives on complexity in the earth system[END_REF][START_REF] Bandt | Permutation entropy: a natural complexity measure for time serie[END_REF]. Some of these quantifiers have been applied to climate science to analyze regime shifts and tipping points [START_REF] Donges | Nonlinear detection of paleoclimate-variability transitions possibly related to human evolution[END_REF]. These methods usually require embedding the time series in a higher dimensional space, a procedure known as "state space reconstruction" [START_REF] Takens | Detecting strange attractors in turbulence[END_REF]. Among these methods the entropy quantifier proposed by Corso et al. [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF] has clear advantages when dealing with large spatiotemporal fields (∼ 10 4 time series in this work given the focus on two-dimensional spatial fields that evolve in time): (a) it avoids the state space reconstruction that may lead to false results [START_REF] Marwan | How to avoid potential pitfalls in recurrence plot based data analysis[END_REF][START_REF] Wessel | Practical considerations of permutation entropy[END_REF] and (b) it relies on just one parameter, the vicinity threshold . Additionally, Corso et al. [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF] showed that given a generic time series, results are (a) well correlated with its maximum Lyapunov exponent and (b) have a weak dependence on the distance threshold .

Here, for a given climate field X(t), we compute the information entropy and its variability for all the time series embedded in a spatial grid to derive its spatiotemporal entropy field S X (t). In order to reduce its dimensionality, we then identify sets of grid cells (domains) that behave homogeneously in terms of their complexity using δ-MAPS [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF], briefly described below.

δ-MAPS

Dimensionality reduction

Given a generic spatiotemporal field, δ-MAPS aims at identifying spatiotemporal patterns, defined as spatially contiguous sets of grid cells with highly correlated activity and referred to as domains. The underlying hypothesis of the domain identification algorithm is that domains have epicenters or cores where their local homogeneity is maximum. The first step is then to identify these cores and expand them. Formally, each grid cell i of a spatiotemporal field S(t), is associated with a time series s i (t). Linear trends are removed from the time series. The similarity between two grid cells i and j is quantified by their Pearson correlation r i,j at lag τ = 0. Cores are identified by defining a K-neighborhood Γ K (i) including a grid cell i and its K closest (in terms of geodesic distance) cells. We then compute the local homogeneity of a cell i, defined as the average pairwise correlation in Γ K (i). A grid cell i is marked as a core if its local homogeneity is a local maximum and greater than a threshold δ. Finally, cores are then iteratively expanded and merged to identify domains [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF]. The total number of domains D will depend on the threshold δ and its minimum sized on the K parameter. Details on the methodology are discussed in Appendix A and more in depth in [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF] and [START_REF] Falasca | Dimensionality reduction and network inference using δ-MAPS: Application to the cesm large ensemble sea surface temperature[END_REF].

Network Inference

The temporal variability of domains is quantified by computing their signals. Given a spatiotemporal field X(t), the signal X A (t) of a domain A is the weighted cumulative anomaly of all time-series within that domain:

X A (t) = |A| i=1 x i (t) cos φ i , (3)
where x i (t) is a time series of length T associated to grid cell i with latitude φ i and |A| is the cardinality of the set. Given D domains, the network is inferred by considering each possible pair of domains A and B and computing their Pearson correlation r A,B (τ ) for a lag range τ ∈ [-τ max , τ max ]. The statistical significance is then tested for the 1 2 D(D -1)(2τ max + 1) correlations adopting the False Discovery Rate (FDR) proposed by Benjamini and Hochberg [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] against multiple testing, and excluding autocorrelations by using the Bartlett's formula [START_REF] Box | Timeseries analysis: forecasting and control[END_REF]. Two domains A and B are connected if there exists at least one significant correlation between the two at any lag in the range τ ∈ [-τ max , τ max ], denoted as R A,B (τ ). If R A,B (τ ) includes the lag τ = 0, the link is left undirected, and treated here as a two-way link A ↔ B. If R A,B (τ ) is strictly positive (negative) the link will be directed from A to B (B to A). Finally a weight w A,B is assigned to each link and based on the covariance between the two signals X A (t) and X B (t) at the lag τ * at which their significant correlation r A,B (τ ) is maximized. In this work, we quantify the degree of connectivity of the climate network, along with its time variability, with the network density metric. Given a weighted, direct network N between D domains, with weights w A,B , we define its density ρ N as:

ρ N = D A,B,A =B |w A,B | D A,B,A =B σ A σ B , (4)
where σ A denotes the standard deviation of the signal X A (t) of domain A. The denominator is never zero and accounts for the case in which all possible connections exist with correlation |r A,B (τ * )| = 1. Therefore, if ρ N = 1 the network is a complete direct weighted graph, with each node connected to every other node and with their weights maximized. If ρ N = 0 the network has no links.

3 Strategy to study climate variability in long transient simulations

Model simulations

We analyze two 6000-years long transient simulations of the global climate system realized with the IPSL model, Vlr01 and Sr02 hereafter. The Vlr01 simulation was run with the standard IPSL-CM5A-LR version [START_REF] Dufresne | Climate change projections using the IPSL-CM5 earth system model: From CMIP3 to CMIP5[END_REF] utilized also in the Coupled Model Intercomparison Project Phase 5 (CMIP5). The atmospheric model resolution is 3.75 • in longitude, 1.875 • in latitude and 39 vertical levels. The ocean model resolution is 2 • , with spatial refinement at the equator and in the Arctic, and 31 vertical levels. The Sr02 simulation has a higher horizontal resolution in the atmosphere (2.5 • in longitude, 1.125 • in latitude) and includes the new 11-layer ORCHIDEE land surface model and a dynamical vegetation module [START_REF] Braconnot | Strengths and challenges for transient mid-to late holocene simulations with dynamical vegetation[END_REF]. The latter allows both the Leaf Area Index and the vegetation type to vary in each grid cell as a function of climate and CO 2 . 1000-years long simulations were performed to produce initial conditions in equilibrium with the external forcing [START_REF] Braconnot | Impact of multiscale variability on last 6000 years indian and west african monsoon rain[END_REF]. Both simulations are forced by yearly-updated Earth's orbital parameters [START_REF] Berger | Long-term variations of caloric solar radiation resulting from the earth's orbital elements[END_REF], and atmospheric composition derived from ice-core reconstructions [START_REF] Braconnot | Strengths and challenges for transient mid-to late holocene simulations with dynamical vegetation[END_REF][START_REF] Otto-Bliesner | The PMIP4 contribution to CMIP6 -part 2: Two interglacials, scientific objective and experimental design for holocene and last interglacial simulations[END_REF]. When showing time-series we always report the simulation time on the x-axis (i.e., time going from 1 to 6000 with 6000 = year 1950).

The focus of the analysis is on sea surface temperature (SST) fields, saved as monthly averages.

Strategy

Given a spatiotemporal climate field, the objectives of our approach are (a) to identify spatially contiguous regions undergoing regime shifts, and (b) to investigate the time dependent connectivity patterns between those regions.

This can be achieved as follows:

• We consider a climate field X(t) and assign to each grid cell i a (constant) value of the vicinity threshold i . The selection of i for a time series x i (t) is based on the rationale that its entropy s i ( i ) should be invariant under small perturbations around its value such that

s i ( i ) ∼ s i ( i ± δ i ), δ i << i .
Details on the heuristics can be found in Appendix B. For every grid cell i we have a (linearly) detrended, anomaly time series x i (t) with T data points. We consider W years of data, and compute its recurrence entropy s i (see section 2.1). We do so every ∆W years and for each i. This step results in the computation of a spatiotemporal complexity field S X (t).

• Given S X (t), its dimensionality is reduced by identifying its domains (see section 2.2.1) and for each domain A the average entropy (complexity) signal is computed as

S A (t) = |A| i=1 s i (t) |A|
, where s i (t) is the entropy time series at grid cell i and |A| the number of grid cells inside the domain. The information on climate regime shifts of each domain A is then contained in S A (t). Changes in information entropy are directly linked to structural changes in the RP and therefore to transitions in the system's dynamics. The domains identified can then be regarded as spatially contiguous sets of grid cells undergoing regime transitions or shifts in a homogeneous way.

• Complexity is a metric that quantifies the variability of a given time series and the D domains can be thought as the true, semi-autonomous components of the system in the 6000 years analyzed (all the grid points within a domain are characterized by highly correlated variability). This allows to consider them as "static" nodes of the climate network and investigate their connections in time, given that the D domains identified can be connected through oceanic and/or atmospheric teleconnections. Their connections are captured by the network topology in the X(t) field. We consider now Y years of data and infer a weighted and direct network between domains. We do so every ∆Y years. This step results in the computation of a time dependent network between all D domains.

This strategy is applied to the two SST fields presented in section 3.1. The time resolution is monthly, and the length of the simulations is 6000 years.

It should be noted that the domain identification and network inference act on different fields: first, the dimensionality of the entropy field S X (t) is reduced by identifying its underlying spatial components, then the connectivity patterns between these domains identified in the entropy field are inferred in the original climate field X(t).

Parameters

We compute the complexity of (linearly) detrended SST anomalies using time windows of W = 100 years every ∆W = 20 years.This choice of W allows for the investigation of periods with decadal and multidecadal variability and shifts between them. We choose ∆W = 20 years as a compromise between having a reasonable time resolution in the entropy field and restricting the large amount of computations. The recurrence entropy of each time series is computed by sampling M = 10000 microstates of size M = 4 (as in Corso et al. [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF]). The dimensionality reduction in the spatiotemporal entropy field S X (t), allows for the identification of regions undergoing the same dynamical transitions. For this step we use a value of α = 0.01 to select the threshold δ (as done in Fountalis et al. [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF] and Falasca et al. [START_REF] Falasca | Dimensionality reduction and network inference using δ-MAPS: Application to the cesm large ensemble sea surface temperature[END_REF]) and K = 2 and K = 6 for the low-and high-resolution simulations, respectively. Results are qualitative similar for other values of K close to the selected ones. Networks between the D domains are inferred using Y = 100 years of (linearly) detrended anomalies every ∆Y = 3 years. This results in a total of 1967 networks computed in each simulation. We choose Y = 100 years for consistency with the entropy computation. A choice of ∆Y = 3 years is feasible given the reduced dimensionality of the dataset (i.e., D domains instead of all grid cells) and has the advantage of quantifying temporal shifts in connectivity patterns with high precision. For the network inference we set τ max = 12 months and and the FDR parameter q = 0.04. Therefore, for every network at most 4% of the inferred links may be false positives.

Results

Mean complexity pattern

We first focus on the mean state of the complexity field in the two simulations. Figure 1a,b shows the mean spatiotemporal entropy field S SST (t) .

The spatial distribution and magnitude of the SST mean complexity is nearly identical across simulations (Figure 1a,b). There is a marked latitudinal gradient between the tropics and mid-to high latitudes (Figure 1c,d) resulting mostly from the large complexity in the equatorial Atlantic and Indian Ocean basins. In the inter-tropical belt, the SST mean complexity depicts significant basin-dependency in both magnitude and spatial gradients. The tropical Indian Ocean is characterized by very large SST mean complexity, with values ranging from ∼ 0.9 (1 = maximum complexity) in its eastern equatorial part to ∼ 0.8 in its western tropical part. Similar values are found in the equatorial Atlantic, also with a small east-to-west gradient. The tropical Pacific mean complexity, on the other hand, exhibits non-trivial spatial features. There is a steep west-to-east gradient in complexity, with the highest values in the upwelling cold tongue and significantly lower mean complexity in its western part. Additionally, the lowest global mean complexity is located in the central Pacific on two hemispheric-symmetric stripes. The equatorial Pacific hosts the strongest mode of interannual climate variability, ENSO, and the potential predictability of tropical SST anomalies is, on average, greater in the Pacific than in the two remaining basins due to the ENSO quasiperiodicity. The analysis suggests that in the 6000 years considered, ENSO-related predictability in the IPSL model is on average larger in the central Pacific, east of 180 • E, in two symmetric bands north and south of the Equator.

In the extra-tropics, the SST variability displays relatively low mean complexity outside the western boundary current pathways and the eastern boundary upwelling systems. The strongest shifts in complexity in both simulations, however, are found at high latitudes, as we will show next.

Exploring climate regime shifts in SST

To investigate the variability in the information entropy field S SST (t), we reduce its dimensionality using δ-MAPS. The identified domains are shown in Figure 2 for both simulations and are similar in terms of their spatial characteristics and mean complexity. Following from Figure 1, domains with large mean complexity are found in the Indian Ocean, Equatorial Atlantic basins, the Southern branch of the Horse Shoe pattern extending from the South Pacific Convergence Zone through Indonesia, and along the path of major boundary currents. We report here on two domains that undergo major shifts in local dynamics in both simulations, to exemplify how the proposed strategy allows to quantify regime transitions in terms of timing and spatial characteristic. The two domains have been chosen because they exhibit considerable jumps in entropy (∆S ∼ 0.1) in periods as short as 200 to 300 years. For each domain we show the time series of SST anomalies in centuries of highest and lowest complexity. A more comprehensive investigation of all regime shifts, their interrelationships and causes is beyond the scope of this work, that focuses on the methodology, and is left for future studies.

First, we focus on the domain in the Vlr01 simulation with the largest regime shift (Figure 3). It connects the Indian and South Atlantic Oceans (Figure 3a) and embeds the Agulhas retroflection current. Its mean complexity is 0.68 (Figure 3b), a value similar to that in the higher resolution simulation in which, however, shifts occur at different times (Figure 3 andAppendix C). Therefore, for this domain, we can conclude that regime shifts are likely due to the internal variability of the system rather than variations in the orbital parameters or gas concentrations. Every point in the complexity signals is computed based on the detrended SST anomalies over 100 years. Panels c and d in Figure 3 show the monthly SST anomalies relative to the century with the lowest and largest complexity, respectively. In roughly 400 years, from year 5120 to 5520, the domain undergoes a sudden, rapid transition from high to low and back to high-frequency SST variability. Specifically, from year 5281 to 5380, we observe the least complex SST variability with a clear multidecadal oscillation (Figure 3c). This multidecadal oscillation has a periodicity of 33.3 years, which is significant at the 95% level according to a red noise test [START_REF] Imkeller | Stochastic climate models[END_REF][START_REF] Dijkstra | Networks in climate[END_REF]. The system then jumps to its highest complexity state in the 6000-year record between year 5421 and 5520 (Figure 3d).

Figure 4 presents the complexity signal for a domain crossing the Drake Passage in the Sr02 simulation. The mean complexity of this region is 0.73, higher than the majority of domains identified south of 30 • S (see Figure 2b). Again, the SST anomalies relative to the centuries with lowest and largest complexity are also shown (Figure 4c,d). Low-frequency oscillations in SST anomalies are observed from year 3361 to 3460 when the lowest complexity state and significant, energetic spectral peaks of 20 years are found. The century between simulation year 1981 and 2080, on the other hand, displays high-frequency variability and corresponds to the highest complexity state. For this domain, the largest jump in complexity appears over a period of about 300 years, from 3200 to 3500 (Figure 4b).

We note that results are constrained by the choices of time windows for the entropy computation (i.e., W = 100 years windows every ∆W = 20 years) and that new insights could be obtained using shorter sliding windows (i.e., ∆W << 20) or longer periods (i.e. W > 100 years). However, regime transitions are also reflected in changes in the topology of the networks between domains, and here, climate networks are inferred every 3 years, allowing for the investigation of regime shifts at very high temporal resolution (see next section).

Exploring variability in SST climate networks

In our approach, the domains in Figure 2 emerge as the true semi-autonomous components of monthly SST variability. They are considered as true components since they are spatiotemporal patterns of the system under study and largely independent of the spatial resolution of the grid and as semi-autonomous since they can interact via coupled atmospheric-ocean dynamics. x-axis: time in years; 6000 = year 1950; the year reported indicates the first year of the period considered i.e., year = 1 for the period going from year 1 to 100).

Here we investigate their time-dependent connectivity patterns using the network density metric ρ N proposed in Eq. 4. Changes in the network density allow for the quantification of transitions in the degree of connectivity on a global scale. We focus on abrupt regime transitions.

We draw attention to two points: (a) a system with no regime transitions would imply constant values of the network density over time and (b) high correlations between the two densities would be obtained if the external forcings (i.e., orbital parameters and trace gases) were the dominant mechanisms for global regime transitions In the 6000 years explored here, the climate system undergoes relatively large dynamical changes in its network density from minima of 0.10 (0.07) to maxima of 0.19 (0.18) for the low (high) resolution simulations, respectively, and the correlation between the two densities is not significant (c.c. = 0.05) (Figure 5). Therefore, both simulations display regime transitions driven mostly by internal variability. The mean network densities are very similar across simulations (i.e., 0.14 and 0.13 respectively for Vlr01 and Sr02). The variability in density, quantified by its standard deviation σ, on the other hand, is greater in the higher resolution Abrupt changes in the variability of the network density are, on average, smaller in the Vlr01 than in the Sr02 simulation, with the maximum from year 1171 to 1270. The density in the Sr02 undergoes frequent, abrupt changes approximatively from year ∼ 2000 to ∼ 4000 and exhibit fewer transitions in the first and last 2000 years of the simulation. The curves are not correlated (correlation r = 0.04), and major transitions in the two climate networks can be considered as an emergent [START_REF] Dijkstra | Numerical bifurcation methods applied to climate models: analysis beyond simulation[END_REF], unpredictable phenomena, arising from the internal variability of the system rather than from the (common) external forcing. There are some exceptions to this view, such as the shifts identified approximately around the years ∼ 1800, ∼ 3200 and ∼ 5800. Around these periods, both simulations display large transitions in global connectivity. Nonetheless, the shifts around year ∼ 3200 is of a different type in the two simulations: around this year, Vlr01 exhibits a large, abrupt increase in connectivity, while Sr02 the opposite (Figure 5). Interestingly, around the year ∼ 1800 (∼ 4200 BP) both simulations undergo a rapid, strong increase in connectivity of the network. The timing of this event agrees with the transition between the so-called Northgrippian and Meghalayan ages in the Holocene [START_REF] Walker | Formal ratification of the subdivision of the holocene series/epoch (quaternary system/period): two new global boundary stratotype sections and points (gssps) and three new stages/subseries[END_REF], that was characterized by severe aridification and cooling that extended from North America to the Mediterranean basin and North Africa, all the way to India and China [START_REF]Cultural responses to climate change during the late holocene[END_REF][START_REF] Staubwasser | Climate change at the 4.2 ka BP termination of the indus valley civilization and holocene south asian monsoon variability[END_REF][START_REF] Rousseau | The 4.2 ka BP climatic event [special issue[END_REF]. This suggests a possible dominant role of the external forcing in causing this specific transition. However, more ensemble members are needed to assess the significance of this result. Furthermore, the investigation of multiple fields including air-temperature and precipitation may shed light on the evolution, connectivity and wide-spread impacts of this regime shift. This is left for future work.

SST complexity in the "real" world

The model simulations analyzed in this study allow to explore SST variability at spatial and temporal resolution unavailable in the observational records. However, climate models are far from perfect, and results are compromised by structural model errors [START_REF] Frigg | Laplace's demon and the adventures of his apprentices[END_REF]. Conclusions may be more informative when presented together with model biases [START_REF] Thompson | Escape from model-land[END_REF] and, therefore, we present a brief, qualitative comparison of the global SST complexity in the IPSL model and in an observational data set. For a first assessment of SST climate networks in the IPSL model we refer the reader to Fountalis et al. [START_REF] Fountalis | ENSO in CMIP5 simulations: network connectivity from the recent past to the twenty-third century[END_REF]. We consider the last 50 years of simulation, from 1900 to 1950, and compare the SST complexity reproduced by the IPSL model in the higher resolution configuration (Sr02) with the COBEv2 reanalysis [START_REF] Hirahara | Centennial-scale sea surface temperature analysis and its uncertainty[END_REF] saved at 1 • degree spatial resolution and as monthly averages (Figure 7). COBEv2 is also a model product and data availability in the first 50 years of the 20th century is rather poor [START_REF] Hirahara | Centennial-scale sea surface temperature analysis and its uncertainty[END_REF]. Nonetheless, we expect the SST variability in COBEv2 to be more realistic than IPSL.

The complexity field shown in Figure 7 for the Sr02 simulation is qualitatively similar to the mean complexity fields shown in Figure 1. This highlights that SST complexity is quasi-stationary and trends are small in the mid-to late Holocene transient simulations. The average latitudinal dependence of complexity in the observational data (Figure 7b,d) suggests larger predictability in the tropics (from ∼ 25 • S to ∼ 25 • N) and lower at higher latitudes (from ∼ 60 • S to ∼ 25 • S and ∼ 25 • N to ∼ 60 • N).This is due to the low entropy values in the equatorial Atlantic and Pacific basins. Therefore, the complexity of the simulated SST in the tropical Atlantic is largely overestimated by the model. Large complexity characterizes the Indian Ocean in both data sets, but is less uniform and weaker in COBEv2 than in the model. The complexity of the ENSO region emerges in the reanalysis as confined to the central and eastern Pacific and is more uniform than in the model. Interestingly, two hemispheric-symmetric stripes marked with low complexity can be identified also in COBEv2 for the ENSO area, suggesting larger predictability potential for the El Niño events just South of the equator and East of 180 • E.

Discussion

This work presents a new framework stemming from dynamical system theory and complex network analysis to investigate regime transitions in spatiotemporal climate fields in a comprehensive way. The approach couples a complexity quantifier [START_REF] Corso | Quantifying entropy using recurrence matrix microstates[END_REF] to a dimensionality reduction methodology [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF] in order to (a) identify spatial regions (i.e., domains) of grid cells with highly correlated evolution of their complexity and (b) study the time-dependent, weighted, direct network between these regions. The approach is applied to two 6000-years long simulations of the IPSL Earth system model. Both simulations develop highly non-trivial spatiotemporal variability, which may stem at time from common forcing (changes in Earth's orbit and trace gases) but more often from intrinsic, internal variability. We tested the proposed strategy focusing on decadal to multidecadal SST variability in the last 6000 years, providing an analysis of regime transitions and time-dependent networks that simultaneously embraces spatial and temporal scales.

The results of our analysis can be summarized as follows:

• SST variability in the IPSL model is characterized by relatively low mean complexity (i.e., high predictability) in high latitudes and high complexity (i.e., low predictability) in the equatorial Indian and Atlantic Oceans, and along major eastern and western boundary currents. The ENSO region displays a non-trivial spatial configuration of mean complexity, with high values in the cold tongue and low values and higher predictability in two hemispheric-symmetric bands in the central Pacific.

• The dimensionality reduction of the entropy field allows for the identification of its underlying domains. Each domain is a set of points with highly correlated variability in information entropy. A domain's complexity signal, in turn, allows for exploring regime shifts in time. We provided evidence for transitions in local system dynamics in the two simulations that are identified as abrupt jumps of large amplitude from low to high-frequency SST variability and vice versa in periods as short as 200 to 300 years. Our results suggest that, at least in the modeled climate, multidecadal oscillations can rapidly emerge and fade due to the system's internal dynamics (i.e., independent of the external forcings). More theoretical and data driven work is needed to better understand the physical underlying mechanisms behind these shifts and their expression and impacts on the overall system and not limited to SSTs.

• The identified domains can be considered as the true semi-autonomous components of the system in the 6000-years long simulation. Domains can interact via coupled atmosphericocean dynamics, giving rise to complex connectivity patterns. These connections, and their variability, can be studied as a time-dependent, weighted and direct network between domains. The time evolution of the network density reveals that the whole system underwent large shifts in the degree of connectivity during the last 6000 years. These shifts occur, in most cases, at different times in the two simulations, implying that most transitions in connectivity may be unpredictable, emergent phenomena [START_REF] Dijkstra | Numerical bifurcation methods applied to climate models: analysis beyond simulation[END_REF]. This suggests that changes in connectivity patterns in timescales as long as 1000 years are dominated by the internal variability of the system rather than by external forcings, at least in the last 6000 years. Finally, we identified a common (among the two simulations) large shift around year 1800 (i.e., 4200 years BP) consistent with previous findings in paleoclimate proxies [START_REF] Rousseau | The 4.2 ka BP climatic event [special issue[END_REF]. A larger ensemble is required to confirm if this shift in connectivity is a response to external forcing. Exploring large regime shifts in different variables and quantifying their underlying causal associations using state-of-the-art methodologies [START_REF] Runge | Detecting and quantifying causal association in large nonlinear time series datasets[END_REF] is left for future work.

• In comparing the SST complexity of the IPSL model with that of a reanalysis product over the first half of the 20th century, differences were mainly found in the tropical Atlantic and Pacific basins. The modeled complexity in the Atlantic around the Equator and in the cold tongue is largely overestimated, while the ENSO domain extends too far west covering most of the equatorial Pacific. We note that evaluating a climate model's ability in reproducing the observed complexity is a very stringent (but useful) test: high similarity could only be obtained if models and observations shared the same underlying dynamics. Such evaluation allows for testing directly the dynamics of the model rather than comparing signal amplitudes or mean states and opens interesting research avenues. Additionally, while dynamical transitions in SST emerge at relatively long-time scales (i.e., hundreds of years), future work could focus on variables expected to show regime shifts at shorter timespans. This would include precipitation or cloud fraction fields in high temporal resolution datasets and/or model integrations.

Transitions or shifts in local and, more importantly, global system dynamics are expected in complex, nonlinear spatiotemporal climate fields. While climate models are more often used for future projections and projected climate shifts cannot be validated, long simulations covering past millennia allow for exploring profound changes in the system dynamics that may have already taken place and for which we have sparse records. In this work we investigated the relationships between these transitions, external forcing and internal climate variability. Our results suggest that shifts in local dynamics and global connectivity patterns in sea surface temperature may often appear as emergent and abrupt phenomena, independently of the external forcing fields because of the internal chaotic dynamics of the system.

A δ-MAPS

Core identification. The domain identification algorithm starts from epicenters or cores. Given a K-neighborhood Γ K (i) of a cell i, we compute its local homogeneity as:

r K (i) = m =n∈Γ K (i) r m,n K(K + 1) . (A.1)
A similar notation is used for the homogeneity of a set of grid cells:

r(A) = m =n∈A r m,n |A|(|A| -1) , (A.2)
|A| being the cardinality of the set and and r m,n the correlation (at lag τ = 0) between the time series embedded in grid cells m and n. A grid cell is marked as a core if its local homogeneity is a local maximum and greater than a threshold δ. More formally, a grid cell i is a core if r K (i) > δ and r K (i) > r K (j), i = j, ∀j ∈ Γ K (i). Cores are then iteratively expanded and merged to identify domains [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF].

Domain identification. Given one or more identified cores c, a spatial grid G and a threshold δ, a domain A is the maximal set of grid cells satisfying three constraint: (a) c ∈ A, (b) I G (A) = 1 and (c) r(A) > δ. I G (A) = 1 denotes that the set A is spatially contiguous. Fountalis et. al [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF] proved that this problem is NP-Hard and relied on a heuristic to solve it. The heuristic algorithm iteratively expands and merges cores to find the full extent of domains. The process starts from the domain with largest homogeneity, let it be A. The expansion algorithm considers all adjacent grid cells of A, and then adds the grid cell i which maximizes the homogeneity r(A ∪ i) and for which r(A ∪ i) > δ. If two domains A and B are adjacent, the merging algorithm determines that the two domains should be merged whenever I G (A ∪ B) = 1 and r(A ∪ B) > δ. This process stops when no more merging and expansions are possible. Two heuristics to infer the δ threshold and the K parameter has been proposed in [START_REF] Fountalis | δ-MAPS from spatiotemporal data to a weighted and lagged network between functional domain[END_REF] and [START_REF] Falasca | Dimensionality reduction and network inference using δ-MAPS: Application to the cesm large ensemble sea surface temperature[END_REF]. Given a significance level α, the threshold δ is computed as a sample average of the statistically significant cross-correlations between randomly chosen grid cells in the dataset considered.

B Heuristic for the vicinity threshold

A crucial step in the computation of a recurrence plot (RP) is the definition of the vicinity threshold . The parameter defines the neighborhood of each state x i in a trajectory of a dynamical system. It must be carefully chosen since a value that is too small may include noisy fluctuations in the RP, and one that is too large could hide the recurrence structure of the time series [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF]. Several "rules of thumb" for the choice of the threshold have been proposed: ∼ 5% of the maximal state space diameter [START_REF] Mindlin | Topological analysis and synthesis of chaotic time series[END_REF], no more than ∼ 10% of the mean (or maximum) state space diameter [START_REF] Koebbe | se of recurrence plots in the analysis of time-series dats[END_REF][START_REF] Zbilut | Embedding and delays as derived from quantification of recurrence plots[END_REF] or a value ensuring a recurrence point density of ∼ 1% [START_REF] Zbilut | Recurrence quantification based lyapunov exponents for monitoring divergence in experimental data[END_REF]. Moreover, choosing strongly depends on the system under study [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF].

Here we present a new simple heuristic tailored for the selection of a different i for each (notembedded) time series x i (t) of a spatiotemporal field X(t).

The proposed strategy relies on the following steps:

1. we consider a monthly spatiotemporal climate field X(t) embedded in a two-dimensional grid at each time step t. For every grid cell i we have a linearly detrended, anomaly time series x i (t) with T data points. We define an i as a percentage ρ of the standard deviation σ i of x i (t). This step defines a new spatial field E(ρ), referred to as -map.

2. For a given -map E(ρ), we compute the RP and its information entropy for every x i (t) in X(t) to get a spatial entropy field S = S(ρ). Since the entropy field depends on ρ, we assess its sensitivity to ρ as follows: first, we compute S(ρ) for ρ ranging from ρ min = 0.01 to ρ max = 1.50 every ∆ρ = 0.01; second, we define a pairwise distance matrix between all S(ρ). Given two fields, S(ρ 1 ) and S(ρ 2 ), we compute their distance as D(S(ρ 1 ), S(ρ 2 )) = i |s i (ρ 1 ) -s i (ρ 2 )|, where s i (ρ) denotes the entropy at grid cell i.

3. Finally, we select an optimal ρ = ρ * such that the computed entropy field does not depend on perturbations around ρ * . This assures that, for every s i ( i ) at grid cell i, s i ( i ) ∼ s i ( i ± δ i ), δ i << i .

In case of results presented in section 5, we obtain a value of ρ * = 0.75 both for COBEv2 and for the Sr02 simulation. We then assign to each grid cell i a value i = ρ * • σ i , with ρ * = 0.75.

For results in section 4, we face an additional problem: we have to compute a spatiotemporal entropy field and not just one time-instant. Therefore, for this step we compute E(ρ) using the complete simulation length. We consider a time series xi (t) of monthly SST defined at grid cell i (with length T = 72000 months) and compute a new time series x i (t) representing the evolution of detrended anomalies at all times t. SST anomalies x i (t) are computed by removing both the seasonal cycle and a linear trend. This is done for each 100-yr non-overlapping window instead of the 6000-yr timeseries to account, at first order, for changes in SST seasonality and (small) trends induced by orbital changes from mid-to late Holocene [START_REF] Braconnot | Impact of multiscale variability on last 6000 years indian and west african monsoon rain[END_REF]. The parameter i is thus constant in time. Given a x i (t), this implies a notion of neighborhood as a function of all possible states explored in the 6000-yr long simulation (i.e., the closeness between two states is conditioned on distances between all states in the 6000-yr simulation). The heuristic is then tested for the first W = 100 years of the simulation. This analysis is found to be independent of resolution. A global minimum in the pairwise distance matrix is identified around

Figure 1 :

 1 Figure 1: (a,b) Mean complexity field for Vlr01 and Sr02 simulations. (c,d) Average latitudinal dependence of complexity in both simulations.

Figure 2 :

 2 Figure 2: (a, b) Domains in the SST entropy field for the lower (Vlr01) and higher (Sr02) resolution simulations, respectively. Colors quantify the averaged complexity values inside each domain.

tFigure 3 :

 3 Figure 3: (a) Domain with the largest regime shift in the low resolution (Vlr01) simulation. (b) Mean complexity signal for the domain in (a) over time (x-axis: time in years; 6000 = year 1950; the year reported indicates the first year of the period considered i.e., year = 1 for the period going from year 1 to 100). (c, d) Associated monthly SST anomalies corresponding to centuries of minimum and maximum complexity, respectively; x-axis: time in years. The black lines indicate a 10-years moving average.

tFigure 4 :Figure 5 :

 45 Figure 4: (a) a domain in the high resolution (Sr02) simulation with a large regime shift. (b) Mean complexity signal over time for the domain shown in (a) (x-axis: time in years; 6000 = year 1950; the year reported indicates the first year of the period considered i.e., year = 1 for the period going from year 1 to 100). (c, d): Monthly SST anomalies corresponding to centuries of minimum and maximum complexity, respectively; x-axis: time in years. The black lines indicate a 10-years moving average.

Figure 6 :

 6 Figure 6: Network density variance in moving windows of 200 years every 3 years (minimum time step) respectively in the Vlr01 (red) and Sr02 (blue). x-axis: time in years; 6000 = year 1950; the year reported indicates the first year of the period considered i.e., year = 1 for the period going from year 1 to 100). Dashed lines indicate periods with qualitative similar shifts at the same time in both simulations.

Figure 7 :

 7 Figure 7: (a,b) Complexity field for the Sr02 simulation and the COBEv2 reanalysis for the 1900-1950 period. (c,d) Average latitudinal dependence of complexity in both fields.

Figure B. 1 :t

 1 Figure B.1: Pairwise distance matrix between the spatial complexity fields computed for the first 100 years as a function of the percentage of standard deviation of the 6000 years detrended anomalies time series defined at each grid cell. Similar results are obtained for the COBEv2 dataset and for other centuries in the simulations.

For a review of abrupt regional shifts and tipping points in the climate system, the reader is referred to the work by Alley et al.[START_REF] Alley | Abrupt climate change[END_REF] and Lenton et al.[START_REF] Lenton | Tipping elements in the earth's climate system[END_REF].

Acknowledgements

This work is a contribution to the JPI-Belmont PACMEDY project (N • ANR-15-JCLI-0003-01). The two climate simulations were performed thanks to a high-end computing PRACE access grand (TRHOL project) and to GENCI-TGCC HPC resources. F.F. thanks Jonathan Donges and Norbert Marwan for suggestions on relevant research in recurrence plots, recurrence networks and entropy quantification, and Gilberto Corso for clarifications on their entropy quantifier. The work was supported by the Chateaubriand Fellowship grant offered by the Embassy of France in the United States.

Data availability

The climate simulations used in this manuscript are freely available (https://doi.org/10.14768/ 20191028001.1). Python codes for recurrence infomation entropy are available at https://github. com/FabriFalasca/NonLinear_TimeSeries_Analysis. δ-MAPS is available at https://github. com/FabriFalasca/delta-MAPS. Additionally, some Mathematica codes for network inference between domains are also available at https://github.com/FabriFalasca/network_Inference.

C Shifts for a domain in the South Atlantic ocean in the Sr02 simulation

Here we present a domain in the Sr02 simulation similar to the one shown in Figure 3