Aofei Mao

Peixun Fan
email: fanpeixun@gmail.com

Loic Constantin

Nan Li

Xi Huang

Bai Cui

Jean-Francois Silvain

Xinwei Wang

Yong Feng Lu

Forming three-dimensional micro-objects using two-dimensional gradient printing

Keywords: controllable 2D-3D transformation, self-shaping, two-photon polymerization, femtosecond laser micro/nano manufacturing, anisotropic shrinkage, stress control TPP two-photon polymerization fs femtosecond

Controllable transformations from two-dimensional (2D) patterns to threedimensional (3D) geometries independent of materials or external stimuli are being pursued in numerous fields. Here, we present an approach to forming various 3D structures through 2D printing using distributed stress inside a polymer. The key is to establish controlled stress fields by introducing composition and property gradients inside a photocurable polymer by femtosecond laser two-photon polymerization. Structural deformation induced by internal stress is a general bottleneck both in materials processing and 3D printing. In contrast to the significant efforts previously made to reduce stress and deformation, we use them to enable shape transformation to construct various 3D micro-objects through 2D printing with engineered stress fields inside. Multi-mode 2D-to-3D structural transformations, including bending, rolling, coiling, waving, spiraling, and out-of-plane distortions are realized in a shape-and locationspecific fashion. This strategy promises a unique way to fabricate delicate 3D objects not feasible through conventional techniques and to circumvent the intrinsic stepping limitations in direct 3D printing using two-photon polymerization. When combined with the standard 2D patterning techniques such as nanoimprint and photolithography, such a 2D-to-3D transformation approach will lay a foundation for high-throughput and cost-effective production of complex 3D nanostructures.

Introduction

Controllable structural transformation from two-dimensional (2D) patterns to three-dimensional (3D) geometries is being pursued in various fields including microelectromechanical systems [1,2], soft actuators [3], biomimics [4], and biomedical devices [5][6][7]. One traditional strategy is to utilize bi-/multi-layer structures composed of different materials, and the mismatches of their properties, such as thermal coefficient and solubility parameter [8,9]. Structural transformation has also been realized using the unique responses of specific materials (e.g., hydrogels [10], shape memory polymers [11], etc.) to external stimuli (e.g., light [12], temperature [13], chemicals [14], magnetic [15], etc.). Origami and kirigami methods have also been developed for structural transformation through the exquisite design of hinges and foldings [16], making use of surface tensions [17,18], cell traction [19], bucklings [20], etc., as the driving forces. However, general approaches capable of realizing controlled 2D-to-3D (2D-3D) transformation to construct 3D geometries with a single material in a single process have yet to be explored.

Direct laser writing via two-photon polymerization (TPP) has been known for maskless and programmable 3D printing at meso-micro-nano scales. However, TPP has primarily been used to directly print 3D designs with little control over the material properties and the subsequent structural transformation. Recent studies reported the reconfiguration and motion of 3D structures after printing by adding stimuli-sensitive compositions (e.g., hydrogels [21], graphene [22], carboxylic groups [23], gold nanorod [24], and bovine serum albumin [25], etc.) into the polymer resins. Spontaneous structural transformation from in-plane 2D patterns to out-of-plan 3D objects independent of specific material compositions can inspire broad interests in micro/nanofabrication but has seldom been achieved.

Here, we use TPP to tailor the local cross-link gradients inside photocurable polymers with submicron resolutions and, thus, to establish spatially distributed stresses for achieving spontaneous, prespecified 2D-3D transformations (Video S1). Internal stresses usually cause undesired structural deformations and, hence, are considered as a common bottleneck both in materials processing and 3D printing. In contrast to the significant efforts previously made to reduce stresses and deformations, we demonstrated that internal stresses can be tailored to form 3D micro-objects instead of directly printing geometries exactly as designed. The printingcondition-determined degree of conversions (DoC) and shrinkage behaviors of polymers were investigated. This research clarifies the critical role of composition and property gradients inside a single material in enabling 2D-3D transformations. Multi-mode transformations were achieved either individually in different shapes or simultaneously in a single complex 3D structure in different photocurable polymers, revealing the universality of the 2D-3D transformation approach. This strategy opens up a new avenue to fabricate 3D micro-objects with complex geometries, reduces the production time cost, and circumvents the intrinsic stepping limitations in direct 3D printing. The fundamentals explored in this study promise new opportunities for addressing general challenges in the fabrication of 3D micro/nanoscale structures for different applications.

Results and discussion

Conception of TPP-enabled 2D-3D transformation

Figure 1a-f illustrate the processes involved in the TPP-enabled 2D-3D transformation approach. The 2D in-plane patterns unfolded from targeted 3D objects are printed via TPP. In TPP, photoinitiators absorb two photons simultaneously in response to femtosecond (fs) laser exposures and decompose into free radicals. The sp 2 carbons from monomers react with surrounding free radicals to form polymer networks until termination occurs. Subsequently, the unreacted initiators and monomers are replaced and removed by developing agents. As the developing agents dry up, the polymer networks tend to shrink [26]. The key step in the TPP-enabled 2D-3D transformation approach is to introduce internal gradients in polymerization. Correspondingly, the level of shrinkage will be location-specific in the 2D patterns, building up stress fields when drying up and enabling the transformation of the 2D patterns to desired 3D objects [27].

As a nonlinear two-photon process, TPP occurs within a nanoscale ellipsoid around the fs laser focal spot (voxel). This enables flexible and accurate control over the stress fields. Specifically, the DoC from monomers to polymer networks directly determines the cross-link densities and, thus, the shrinkage at different locations [28,29]. Since the DoC is laser-dose dependent, the shrinkages and the consequent stress fields can be tailored voxel-by-voxel through tuning the TPP printing parameters, including the hatching distance (HD), slicing distance (SD), scan speed (SS), and laser power (LP). Therefore, the structural transformation from in-plane 2D patterns to out-of-plan 3D objects can be controlled. More details of the 2D-3D transformation approach and its mechanisms are shown in Figs. S1-S3.

As an example, Fig. 1g shows a 3D standing flower structure composed of a pedicle, multiple petals at three different height positions and several stamens, which was fabricated through the 2D-3D transformation. Different DoC gradients were introduced at different locations during TPP printing. Specifically, the pedicle was printed such that one side was polymerized to a higher degree than the other side. For the stamens, higher and lower DoC regions varied spirally during the printing. In the petals, the DoC gradients were set along both the thickness and the radial directions to achieve different curvatures. We succeeded in simultaneously realizing different types of spontaneous 2D-3D transformations in the pedicle, stamens and petals while keeping the entire flower free-standing. The stamens twisted to show different helical shapes (Fig. 1j). Meanwhile, some petals bent upwards while others bent downwards (Fig. 1h,I,k). The head of the pedicle bent towards the lower left together with the stamens and petals, which can be more clearly seen from the zoomed images in Fig. S4b. The dark-field and fluorescence images demonstrate that the 3D flower is intact and undamaged during the transformation (Fig. 1l,m). Besides, Fig. S4a shows that an array of the 3D standing flower were identically fabricated in one batch, confirming the reproducibility of the 2D-3D transformation.

Engineering stress field for directional 2D-3D transformation

By printing individual cubic blocks with different TPP parameters, the "laser dose-DoCshrinkage" relationships were systematically analyzed. The DoC can be determined by measuring the ratio of the C=C to C=O bonds before and after TPP using Raman spectroscopy [30]. The volume shrinkage (S v ) was estimated by measuring the volume difference between the printed and designed cubes and calculated as:

, ( 1
)
where V 0 is the cube volume designed (i.e., 50×50×50 µm 3 here), and V is the cube volume measured after development. When the HD and SD are fixed, a higher laser dose (e.g., a higher LP or a lower SS) results in a lower shrinkage and a higher DoC (Figs. S5-S7). Therefore, a negative correlation between the shrinkage and DoC is revealed, no matter either the LP or SS was used as the controlling factor (Fig. 2a-b). Besides, larger HDs result in higher shrinkages even that similar DoCs are reached, which becomes more obvious at lower DoC regions.

It was reported that the DoC of the IP-Dip resin (used in this research) can be correlated to the TPP parameters as [31]:

, ( 2
)
where c 1 ρ 0 = 0.62, c 2 = -0.15, δ 2,eff = 3 × 10 -55 cm 4 s, c 0 is a constant representing the relevant optical parameters used in TPP, and D v = c 3 / HD n with c 3 = 1.5e 4 μm -2. 36 and n = 0.64. Using the TPP parameters-DoC and DoC-shrinkage correlations as a parameter-window indicator, the shrinkage of TPP samples can be controlled via controlling DoCs under different TPP printing conditions. Further, the shrinkage differences among different portions of the TPP samples can be designed and used to establish programmed stress fields. Through this effort, the stress fields can be utilized to control 2D-3D transformation rather than a detrimental factor that usually causes undesired deformations in TPP.

For 2D-3D structural transformation, anisotropic shrinkages are favored to gain directional control, which is essential for producing complex 3D structures from 2D patterns [16,[32][33][34]. Our experimental results show that the HD plays a key role in the shrinkage anisotropy (Figs. 2cd, S8 and S9). When small HDs (e.g., 0.1-0.3 µm) are used, the shrinkages are similar in both crossing directions, i.e., parallel ( ) and perpendicular ( ) to the hatching direction. In contrast, using a large HD (e.g., ≥ 0.4 µm) results in obviously higher shrinkages in the direction perpendicular to the hatching lines. Moreover, the shrinkage anisotropy can be tuned by the TPP parameters, e.g., LP and SS. When the HD is above 0.3 µm, the shrinkage anisotropy becomes more obvious with lower laser doses (e.g., lower LP or higher SS).

The anisotropic shrinkage originates from the Gaussian profile of the fs laser used, which results in DoC gradients along the radial direction in the cross-section of each polymer line (Fig. S10). During the overlap of TPP lines to cover a volume, the DoCs in the overlap regions between adjacent lines significantly determine the shrinkage behaviors. Using large HDs, lower DoCs are generated in the overlap regions than in the center of each line, making the overlap regions shrink more (Fig. 2e). Using lower LP or higher SS, polymer lines with smaller crosssections are printed, which also causes lower DoCs in the overlap regions (Fig. 2f). Both conditions can make the whole printed volume inhomogeneous, with more shrinkage in the direction perpendicular to the TPP lines. When multiple layers are printed with different HDs, LPs or SSs, directional stress fields with controlled magnitudes can be established to enable 2D-3D transformation.

As a quantitative demonstration, Fig. 2g shows the SEMs of blocks (50×50×50 µm 3 ) printed with all the parameters fixed (LP = 19.5 Mw, SS = 10 mm/s, HD = 0.6 µm, and SD / HD = 2) other than the hatching direction rotated from horizontal (0º) to vertical (90º). The shrinkages are anisotropic and depend on the hatch direction. Further, taking the bending transformation as an example, the anisotropic shrinkage in thin films (160×160×6 µm 3 ) with DoC gradients intentionally introduced along the thickness direction was investigated (three layers were printed by setting HD/SD at 0.4/0.8, 0.3/0.6, and 0.15/0.3 µm, respectively, in the top, middle, and bottom layers). Theoretically, if the low DoC layers shrink isotropically, the thin films would bend spherically because the second-moment inertia along both directions parallel and perpendicular to the hatching direction are equal due to identical geometries. However, the films rolled up with their axes parallel to the hatching directions (Fig. 2h), which indicates that the shrinkages of the low-DoC layers mainly occur in the direction perpendicular to the hatching direction. Figure 2h also verifies that the bending direction of the polymer films can be precisely and continuously tuned by simply changing the hatching direction. Therefore, the hatching direction can be used to control the path of the 2D-3D structural transformation provided that large HDs introduce sufficient DoC gradients (Fig. S11).

In addition to the bending directions, bending curvatures can also be controlled by tuning the anisotropy in shrinkage and the shrinkage difference between the adjacent TPP layers. To gain a deeper understanding of the 2D-3D transformation, we developed a simplified model to evaluate the bending trend of the TPP 2D films, where the 2D patterns can be equivalent to a series of local pairs of high-and low-DoC layers (Fig. S12a-b). It appears that the DoC differences in two adjacent layers determine the bending curvatures of the TPP films with the same thickness. Larger DoC differences result in sharper bending (Fig. S12c-g). The bending-DoC relation can guide the design of local bending to eventually control the overall 2D-3D transformations.

2D-3D transformations to form diverse 3D structures

With the programmable anisotropic shrinkage and stress fields, we can flexibly engineer the polymer shrinkage by TPP and achieve diverse 3D structures transformed from 2D patterns. In addition to the basic single-direction bending (namely rolling) in Fig. 2h, we have also realized coiling transformation by generating DoC gradients in a 2D parallelogram pattern, forming a microcoil with an even pitch and diameter (Fig. 3a). The hatching lines were made parallel to the short edges of the parallelogram and the HDs increased from 0.2 µm at the bottom to 0.6 µm at the top. In contrast, the same 2D pattern without DoC gradients remained flat after drying. The advantages of free design in TPP can be applied to fabricate most complex structures. For instance, by alternately reversing the DoC gradients at different lateral positions in a 2D pattern, waving transformation was realized, producing a bi-curvature structure, such as a wavy strip (Fig. 3b).

Since the shrinkage is primarily determined by DoC, the DoC gradients can also be used to generate out-of-plane distortions and even freeform transformations, leading to 3D structures with dual-curvatures and ultimately spherical geometries. This is beyond what can be achieved by simply bending and folding polymer films using conventional 2D-3D transformation routines [35,36]. For example, we set the laser scan paths to concentric rings when printing a 2D circular pattern and created DoC gradients along both the thickness and radial directions by changing the HD and SD (Fig. 3c). After drying, the flat pattern was spontaneously transformed into a cap shape. In comparison, the 2D circular pattern remained flat when no gradients were introduced.

To further demonstrate the capability of the 2D-3D transformation, microflowers were also achieved by introducing DoC gradients along both the thickness and radial directions (Fig. 3d). As a more general demonstration, we printed a microfiber with a diameter of 6 µm and a height of 200 µm with the DoCs designed to have a spiral-like distribution (Fig. 3e). The straight fiber was transformed to a helical shape after drying. We further explored the possibility of integrating different types of transformations in a single complex structure and succeeded in producing an entire standing flower(Fig. 1g), verifying the flexibility of the 2D-3D transformation to fabricate objects with structurally complex components.

Besides, by adding specific features on the TPP 2D films, 3D micro-objects having applications in different fields can be fabricated via the 2D-3D transformation. For example, we have achieved a closed 3D microstent transformed from a 2D grid pattern by precisely tuning the DoC gradients along the thickness direction, and a 3D micro-origami structure from a 2D threepointed-star plate by creating DoC gradients in the hinge regions (Fig. S13). The 3D microstent and micro-origami structures have potential applications in biomedicine and microelectromechanics, respectively. Additionally, a 3D curved truss structure was formed through the transformation of a 2D flat pattern composed of multiple hexagon cells (Fig. S14). The spaces within and among the cells provide more tunability on the 2D-3D structural transformation.

Shape-and location-specific 2D-3D transformation

Furthermore, the 2D-3D transformation is quantitively controllable, making shape-and locationspecific transformations possible. The shape-and location-specific transformation can be used to mimic unique morphologies and different growing stages of natural plants. For example, microcoils with two to six turns and diameters of 60-240 µm were formed from 2D patterns of the same size (Fig. 4b-d), mimicking the curls of the albuca namaquensis grass (Fig. 4a). Different bending effects were also realized through gradually increasing the DoC gradients along the radial direction (Fig. 4f-i), mimicking the lotus flowers in different blooming stages (Fig. 4e). The 2D-3D transformation can also tailor a specific one or more petals of the flower structure (Fig. 4j-l) and occur to different degrees at each of the eight petals (Fig. S15). Moreover, microhelices with different spins were fabricated by manipulating the DoC gradients in periodicities of π, 2π, 4π, and 8π at a fixed height (Fig. 4n-q), mimicking the spiraling stems of the juncus effuses grass (Fig. 4m). Microhelices with different spins achieved by manipulating the 3D gradient of DoC periodically.

Advantages of 2D-3D transformation over direct 3D printing

Additionally, 2D-3D transformation can fabricate complicated and delicate 3D objects which easily collapse, are flawed, or require special protocols in direct 3D printing [37]. Figure 5 shows that the 2D-3D transformation achieves microcoils with obviously smoother surfaces and improved integrity, while the directly printed 3D structures show obvious stepping features. The intrinsic layer-by-layer manner in direct 3D printing leads to inevitable stepping features on curved surfaces. The unpredictable instability of over-hanging features during direct 3D printing causes failure in fabrication, imperfect stitching of adjacent parts, and more surface roughness (Fig. S16). Both intrinsic limitations in direct 3D printing are circumvented using the 2D-3D transformation. Furthermore, the 2D-3D transformation fabricates complex objects in shorter times, e.g., microcoils and microstents are fabricated three to four times faster than direct 3D printing (Fig. S17). All results above were achieved using the IP-Dip resin. Also, we have tested implementing our approach with other photoresists, including the IP-S resin and, more generally, the SU-8 resin. Both resins have higher Young's modules than the IP-Dip resin, making the structural transformations more difficult. Nevertheless, we succeeded in using both resins to fabricate the entire standing flowers (Fig. S18), verifying the material independency of our approach. Besides, we believe that the concepts demonstrated in this study can be important for addressing general challenges (e.g., deformation) in the fabrication of 3D micro/nanoscale structures designed for different devices and functions [38]. One representative example is the fabrication of 3D curved lattice structures which are major components of architectured metamaterials and micromachines [39,40]. Usually, such 3D lattice structures can only be fabricated via 3D printing. We further demonstrated the capability to flexibly control their curved orientations and spatial topologies which are challenging in direct 3D printing (Fig. S19). By setting appropriate DoC gradients during printing, the lattice structure transformed to a bridge-shaped geometry with all lattice features well maintained. The capability to control the spatial orientations and topologies of the 3D lattice structures can inspire new insights and designs of novel architectured metamaterials and micromachines for various functionalities. In addition, when combined with advanced beam focusing and scanning systems, such as micro-lens arrays, spatial light modulators, and other multi-foci techniques, the throughput of this approach can be significantly improved to meet the needs in practical applications [41].

Conclusions

In this research, we present a general, controllable, and flexible strategy for 2D-3D structural transformation at the microscale, which lays a foundation for constructing diverse, complex 3D objects from 2D patterns. We succeeded in establishing spatially distributed stresses by tailoring the spatial gradients of cross-links in photocurable polymers via TPP. The relationships among the polymerization conditions, polymer shrinkages, and internal stresses were investigated. Multi-mode 2D-3D transformations, including bending, rolling, coiling, waving, spiraling, and out-of-plane distortions, were realized either individually or simultaneously in a single process. Various 3D geometries, such as micro coils, waves, caps, flowers, helices, and micro-origamis, were demonstrated with their transformation states and geometries being controlled precisely, mimicking unique morphologies and different growing stages of natural plants. Such a 2D-3D transformation approach also circumvents the intrinsic limitations of stepping and stitching in direct 3D printing via TPP, achieving 3D objects with much smoother surfaces and improved integrity.

Materials and methods

Two-photon lithography of 2D layered patterns: Two-dimensional polymer layered patterns were fabricated using a commercial TPP system (Photonic Professional GT from Nanoscribe GmbH) with three commercial photoresists: the IP-Dip and IP-S resins from Nanoscribe GmbH, and the SU-8 resin (SU-8 2005) from MicroChem. The dip-in laser lithography (DiLL) configuration was used for printing with the IP-Dip resin, while the oil immersion configuration was used for printing with both the IP-S and SU-8 resins. The fs laser was focused by a Carl Zeiss Plan-Apochromat objective lens (63× NA1.4 Oil DIC M27) and scanned by a computercontrolled galvo in the resins to prepare arbitrary 2D layered patterns with designed DoC gradients. The specific printing conditions for bulk shrinkage tests with the IP-Dip resin are summarized in the Supplemental Table . For the SU-8 resin, prebaking (sequentially at 50 ºC for 5 min, 65 ºC for 30 min, and 95 ºC for 75 min) was applied, followed by cooling down naturally to room temperature before TPP printing.

Sample development:

Both entire standing-flower structures printed with the IP-Dip and IP-S resins were immersed in propylene glycol monomethyl ether acetate (PGMEA) (3M Company) for 1 h to thoroughly dissolve the unreacted resin after the TPP printing, followed by a rinse in 2propanol (Sigma-Aldrich) for 20 min. The entire standing flower printed with the SU-8 resin was rinsed in PGMEA for 12 h and then in acetone (Sigma-Aldrich) for 20 min. To prevent the collapse of their structures due to capillary forces, the standing-flowers were subsequently dried using a critical-point dryer (Tousimis Autosamdri 931) with slow-fill pressure, fill time, purge time, post purge time, critical point time, and vent pressure set at 700 PSI, 150 s, 10 min, 90 s,

Figure 1 .

 1 Figure 1. Overview of the spontaneous 2D-3D structural transformation. a-f) Schematics of the processes of the 2D-3D transformation approach. g) An overall SEM image of a standing flower after the 2D-3D transformation. h-k) Zoomed-in SEMs of the stamens (j) and petals at three different heights (h, i, k). l) Darkfield optical microscopy, and m), green-fluorescence images of the standing flower after transformation.

Figure 2 .

 2 Figure 2. Controlled anisotropic shrinkage through programmed DoC and directional 2D-3D transformation. a, b) Relationships between volume shrinkages and DoC with LP (a) and SS (b) as the controlling factors. c, d) Evolution of shrinkages with LP (c) and SS (d) at different HDs in both crossing directions. e, f) Schematics of anisotropic shrinkage and the consequent bending transformation perpendicular to the voxel lines through

Figure 3 .

 3 Figure 3. Diverse structures achieved through multi-mode 2D-3D transformations. a) A coil structure formed by creating DoC gradients along the thickness direction. b) A wavy structure formed by creating periodic DoC gradients along the lateral direction. c, d) A cap I and a flower structure (d) formed by creating DoC gradients along both the radial and thickness directions. e) A helical structure formed by manipulating DoC gradients.

Figure 4 .

 4 Figure 4. Shape-and location-specific 2D-3D transformation. a) Albuca namaquensis grass. b-d) Microcoils with different pitches and diameters achieved by varying DoC along the thickness direction. e) Lotus flower. fi) Microflowers at different blooming stages achieved by varying DoC along both the thickness and radial directions. j-l) Targeted transformations at one, two, and four petals. m) Juncus effusus spiralis grass. n-q)

Figure 5 .

 5 Figure 5. Advantages of the 2D-3D transformation approach over direct 3D printing. a, b) and d, e) SEM images of the same microcoil structure fabricated through direct 3D printing and 2D-3D transformation, respectively. b) and e) are the magnified images of a) and d), respectively. c) and f) are profiling results at the surface areas marked by red dashed lines in b) and e), respectively.

Data Availability

The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations. However, most materials and processing parameters and conditions are provided in the Supplementary Materials.

120 s, and 550 PSI, respectively [42]. The other structures printed with the IP-Dip resin were also immersed in PGMEA for 1 h at first before immersed in NOVEC-7100 engineering fluidics (3M Company) for 10 min and left to dry naturally under the ambient condition.

Sample characterization and shrinkage measurement:

A Hitachi S4700 field-emission scanning electron microscope (FE-SEM) was used to observe the shapes and morphologies of the 3D structures after the transformation. The dark-field and green-fluorescence images in Fig. 1 were collected using an Olympus AX70 microscope with either a 20× or 40× objective lens. The excitation wavelength of ~550-570 nm and the emission of fluorescent light in a wavelength range of ~580-600 nm were used for fluorescence imaging [43]. The dimensions of the printed samples were measured by ImageJ software on their SEM images, through which their volume and directional shrinkages were calculated. Doc measurement: Raman spectroscopy was used to determine the DoC of the TPP reaction. Raman spectra were acquired by a Renishaw inVia H 18415 Raman microscope. The excitation source was an argon laser with a wavelength of 785 nm and a power of 50 mW. The excitation and signal collections were carried out using a 50× objective lens with a numerical aperture of 0.50. Raman spectra were recorded using an accumulation time of 15 s. The DoC values were calculated from the Raman peak intensities of C=C and C=O as:

(3) where A C=C , A C=O and A' C=C , A' C=O are the integrated intensities of the corresponding peaks in the TPP structures and the unpolymerized resin, respectively. Besides, a Nikon A1R confocal fluorescence microscope was used for mapping DoC under different laser doses using the Nikon NIS-Elements imaging acquisition and analysis program.

Profile measurement: The structural profiles were measured by a Keyence laser scanning microscope (VK-X200K) which works at a laser wavelength of 408 nm. A Nikon objective lens (CF Plan Apo, 150×, NA 0.95) was used to capture the fine details on the structure surfaces.

Supplementary data

Supplementary data to this article can be found online.

Conflict of Interest

The authors declare no conflict of interest.

Acknowledgements


Léa Laurent
email: lea_laurent02@etu.u-bourgogne.fr

Albin Ullmann
email: albin.ullmann@u-bourgogne.fr

Claude Perrot
email: cperrot@groupama-ra.fr

Xavier Aubouy
email: xavier.aubouy@groupama.com

Thierry Castel
email: thierry.castel@u-bourgogne.fr

How abrupt change in surface temperature impacts water cycle over France ?

Keywords: climate hazard, water balance, abrupt shift, warming aléa climatique, bilan hydrique, rupture abrupte, réchauffement

Along with the rapid warming tendency of the last decades over France, a shift in air temperature has been detected in 1987/88, delimiting two main climate periods from 1959 to 2019. Evidence of pre-and post-shift runoff seasonal differences suggest an impact of this abrupt warming on water cycle evolution. As yearly precipitation slightly evolve during the studied period, evapotranspiration appears to be one of the main drivers of the hydrological cycle evolution. Spring is marked by a strong increase in evaporative demand, followed in summer by a rise in water stress and a drying of soil water content. Besides temporal disparities, this local water cycle evolution depicts different spatial responses over France. Modifications of climate hazard linked to water cycle is expected to impact agro-climatic risks and crop yields.

Introduction

Since 1980s in France, the warming trend intensifies strongly, consistent with climate simulations including the anthropogenic forcing [START_REF] Terray | Quantifying 21st-century France climate change and related uncertainties[END_REF]. As a result of this warming tendency, an abrupt shift is detected [START_REF] Brulebois | Hydrological response to an abrupt shift in surface air temperature over France in 1987/88[END_REF][START_REF] Reid | Global impacts of the 1980s regime shift[END_REF] in maximum and minimum surface air temperature around 1987/1988, delimiting two main climate periods: 1959-1987 and 1988-2009. Comparison of runoff pre-and post-shift shows significant decrease between January and July. Variations in runoff response suggest local water cycle changes with a "soil memory effect" : negative soil moisture anomalies in spring lead to drier soils in summer, inducing lower runoff at the end of summer [START_REF] Boé | Multi-decadal river flow variations in France[END_REF]. Which climatic variables, i.e. precipitation and/or evapotranspiration, are responsible for this local water cycle evolution? Beyond that, what are the spatial and the temporal pattern changes of water balance? As there is no necessary correlation between temperature changes and long term water cycle and water balance variations (Senevitrane, 2012), it is necessary to highlight this subtle balance.

Material and methods

Study area and data

The initial dataset used for this study comes from the Safran-Isba-Modcou model (SIM2) developed by Météo France [START_REF] Soubeyroux | Safran-Isba-Modcou (SIM) : Un outil pour le suivi hydrométéorologique opérationnel et les études[END_REF]. This 8-km-resolution reanalysis chain combines a meteorological analysis system, a land surface scheme and a hydrogeological model. A validation of the whole SIM chain show that the model is quite robust both in space and time, and gives a good estimation of the water fluxes [START_REF] Habets | The SAFRAN-ISBA-MODCOU hydrometeorological model applied over France[END_REF]. The quality of the SIM chain for the representation of the actual climate at a fine resolution is a major asset for climate change impact studies [START_REF] Soubeyroux | Impact du changement climatique sur la sècheresse et l'eau du sol en France : les résultats du projet CLIMSEC[END_REF], Vidal et al., 2010b). With an 8km spatial resolution at daily time step, SIM interpolated precipitation and computed PET (Potential EvapoTranspiration) and AET (Actual EvapoTranspiration) from 1959 to 2019 data are relevant to address the complexity of processes leading to changes in local water cycle [START_REF] Soubeyroux | Safran-Isba-Modcou (SIM) : Un outil pour le suivi hydrométéorologique opérationnel et les études[END_REF]. Based on the surface warming abrupt shift, these data allow us to quantify the evolution of climate hazard linked to local water cycle on a continuous time-scale and over the entire French territory.

Agro-climatic indices

Four variables among all SIM outputs are considered here: maximal 2m-temperature (Tmax), liquid precipitation (PRELIQ), potential and actual evapotranspiration (PET and AET). These daily variables are aggregated to monthly, seasonal and annual scale and are used to compute several types of indices providing informations on the evolution of water balance and hydrological cycle components.

Potential evapotranspiration (PET) [START_REF] Doorenbos | Guidelines for predicting crop water requirements[END_REF], considered as a climatic measurement, is introduced to study the evaporative demand of the atmosphere only affected by climatic parameters [START_REF] Allen | Crop evapotranspiration -guidelines for computing crop water requirements[END_REF]. When cultivating crops in fields, the agronomic management and environmental conditions can differ from standard conditions, leading to an actual evapotranspiration (AET) that may vary from PET. In the SIM chain, PET and AET are computed in Isba land surface scheme using methods and equations described in [START_REF] Noilhan | A Simple Parameterization of Land Surface Processes for Meteorological Models[END_REF] and in [START_REF] Noilhan | The ISBA land surface parameterisation scheme[END_REF]. Isba model as it is used in the SIM chain is driven by incoming radiation, precipitation, atmospheric pressure, air temperature and humidity, and wind speed at a reference level. Along with PET and AET, two indices depicting the evolution of the water balance are computed: the difference PET -AET informs on the soil water deficiency, the ratio AET/PET depicts the efficiency of the evapotranspiration process. The study of the evolution of those two indices allows to assess the impact of temperature increase on the evolution of water balance and crops hydric stress.

Statistical analysis

The four daily series of the variables of interest are aggregated to monthly, seasonal (winter: December-January-February (DJF); spring: March-April-May (MAM); summer: June-July-August (JJA); autumn: September-October-November (SON)) and annual time scale. Robust statistical tests such as change point detection and distribution comparison techniques are used to assess the characteristics of the evolution in precipitation, PET and AET due to temperature evolution.

As change-points detection techniques were used in [START_REF] Brulebois | Hydrological response to an abrupt shift in surface air temperature over France in 1987/88[END_REF], they are used here to identify shifts in climate data series. Among all techniques available to study climatic shifts, Bayesian approaches appear to be adapted to the characterization of abrupt changes [START_REF] Barry | A Bayesian Analysis for Change Point Problems[END_REF], especially in climatic records [START_REF] Ruggieri | A Bayesian approach to detecting change points in climatic records[END_REF]. The method assesses posterior mean and change-point probability, and quantifies the level of confidence when change-points are identified as significant.

Significant differences in mean values for liquid precipitation and PET and AET before and after 1987/1988 temperature shift are detected with a robust Bayesian test [START_REF] Kruschke | Bayesian estimation supersedes the t test[END_REF]. This estimation is based on the Bayesian posterior probability distribution, evaluating whether the probability of a difference is high enough to matter. By using the values of the 95% Highest Density Interval (HDI) to define the confidence interval, the Bayesian method provides information about the magnitude and significance of the difference between two distributions.

Results

Detection of shift in temporal evolution of variables linked to the hydrological cycle

Temporal analysis of variables linked to the hydrological cycle depicts important shifts in PET and PET -AET at the end of the 1980s, following the evolution pattern of maximal temperature (Figure 1). Yearly, the mean of the 1988-2019 period for PET is around 100mm higher than the mean of the 1959-1987 period (720mm) (Figure 1a). AET increases to a lesser extent than PET (+50mm in average between the two periods), resulting to an abrupt increase in PET -AET from 200mm before the 1987 shift up to 250mm after. As a result of a sharp increase in maximum air 2m-temperature in 1987/1988 over France, evaporative demand (PET) increases. Vegetation cover and soils are not able to respond to this rise at the annual scale (AET), revealing the large deficit of soil water content in the second period, implying an increase in hydric constraints (PET -AET). Yearly precipitation remain quite stable during the whole period, confirming the major influence of temperature on the evolution of water demand and soil water content.

Contrasted seasonal changes highlight a strong increase in evaporative demand in spring after 1987/1988 (Figure 1b), followed by a rise in hydric constraints. Our results show that water constraints strongly increase in summer (Figure 1c), along with the drying of soil water content.

In autumn (Figure 1d) and winter (Figure 1e), water demand and hydric constraints increase more slightly or even decrease after the shift, which could be explain by the refill of soil water reservoirs with precipitation.

Since 2010 however, a new step of water constraints seems to be reach as shown with the increase of the PET -AET index values. Summer and even autumn depict the same pattern, pointing out the extension of the drying period of soils during the year.

1987/1988 shift: a turning point

The evolution of the means of the hydrological cycle variables before and after the 1987/1988 shift depicts annual and seasonal temporal disparities (Tables 1, 2 and3). PET and AET means both increased significantly on the 1988-2019 period compared to the 1959-1987 period at the annual scale. This evolution is larger and more widespread for PET (+90mm, 45 % of SIM grid points) than for AET (+45mm, 9 % of stations), highlighting the rise in evaporative demand and in water stress after 1987/1988.

Spring experiences the highest percentage of stations recording a shift in 1986/1987/1988 for PET and AET (Tables 1 and2), directly impacted by the brutal increase in air 2m-temperature in this season [START_REF] Brulebois | Hydrological response to an abrupt shift in surface air temperature over France in 1987/88[END_REF]. Summer display the larger gap between PET and AET anomalies, rising from respectively 39mm and 8mm. Precipitation slightly increases yearly and in summer and autumn in the second period, but always to a lesser extent than PET (Table 3). These results point out an important increase in evaporative demand in spring, draining soil water content after 1987/1988 shift. Evaporative demand still rises in summer, resulting in an aggravation of water constraints between June and August, even extending into autumn. Autumn and winter are characterised by a refill of the reservoirs, although a rise in water constraints is depicted in the recent decade during autumn season. As for precipitation, the evolution of the variables linked to the hydrological cycle outlined in these results may hide geographical disparities.

Geographical response of water balance to temperature shift

Anomalies of the variables and indices linked to the hydrological cycle depict spatial diversity in their evolutions (Figure 2). PET is rising almost everywhere after the 1987/1988 air 2m-temperature shift, confirming the strong influence of the abrupt warming on the evolution of water demand. Vosges, Jura and Morvan massif, centre of France, Ardèche, Cévennes and Haut-Languedoc regional natural reserves and southwest of France experience the strongest increase in evaporative demand, the mean of the 1988-2019 period being around 30% higher than the mean of 1959-1987 period. AET is also almost increasing everywhere, with a 30% rise after the 1987/1988 shift in the Vosges, Jura and Massif Central. Even if some contrasted evolution are show at local scale, the south of France shows roughly a decrease of liquid precipitation whilst north part (above 45N°) shows a clear increase.

Consequently, means of the AET/PET index present contrasted territorial evolutions. On the extreme northeast, in Picardie, along the Rhône valley, on the mediterranean coast and southwest of France, PET increases more than Seasonal evolutions of PET and AET depict various geographical changes (not shown), affecting the evolution of local hydric constraints in different ways throughout the seasons. In spring, AET/PET ratio decrease up to 20% over the entire France except mountainous regions, highlighting the beginning of the water stress increase. These hydric constraints extend on almost all the French territory in summer, AET/PET ratio decreasing up to more than 30% on some areas.

Most of the regions experiencing an important increase in hydric constraints after the 1987/1988 shift in air 2m-temperature at the annual and seasonal scales are large agricultural production areas.

Conclusion

As expected, our results clearly point significant changes in the elements of the local water cycle. One interesting result shows that, at yearly scale, a large part of the territory has been able to follow the large increase in PET. However, subtle evolutions are shown seasonally and at local scale. Summer is particularly impacted by an increase in hydric stress, PET -AET being around 30mm higher after the 1987/1988 shift. Extreme northeast of France, Picardie, Rhône valley, Mediterranean coast and southwest of France are the most exposed territories, AET/PET ratio rising between 10 to 20% on the post-shift period at the annual scale. Furthermore, it appears that the last decade seems to experiment a new acceleration stage of the hydrological cycle that may affect water balance. This latter is of primary importance for the productivity of grass, crop and forest covers [START_REF] Creed | Climate warming causes intensification of the hydrological cycle, resulting in changes to the vernal and autumnal windows in a northern temperate forest[END_REF][START_REF] Sergent | N. Douglas-fir is vulnerable to exceptional and recurrent drought episodes and recovers less well on less fertile sites[END_REF]. Since the 1990s in France, a trend towards stagnation of crop yields is observed, bread wheat being particularly affected [START_REF] Schauberger | Yield trends, variability and stagnation analysis of major crops in France over more than a century[END_REF]. The evolution of climate hazard linked to the water cycle depicted in this study is thus expected to be one of the main causes of crop yields evolution. Hence, consecutive agro-climatic risks evolution is of major concern for decision-makers, and may lead to adaptation process from managers.

Figure 1 .

 1 Figure 1. (a) Yearly and (b, c, d, e) seasonal evolution of variables linked to hydrological cycle between 1959 and 2019. PET: Potential EvapoTranspiration; AET: Actual EvapoTranspiration. Gray bars: SIM chains values; gray solid curves: Bayesian posterior mean values. Solid coloured curves: Bayesian posterior mean values; dashed coloured curves: SIM chains values. Vertical dashed line indicate 1987.

Figure 2 .

 2 Figure 2. Yearly relative anomalies of period 1988-2019 compared to period 1959-1987 for all variables linked to the hydrological cycle.

Table 1 .

 1 Percentage of SIM grid points recording a shift in 1986/1987/1988, anomalies (ΔPET) and 95% confidence interval of the difference in the mean (95% HDI) at annual and seasonal scales, for potential evapotranspiration (PET).

	Period	Potential evapotransiration (PET)
		% of stations	ΔPET (mm)	95 % HDI
	year	45	+90.3*	89.8 -90.9
	DJF	31	+5.73*	5.64 -5.81
	MAM	69	+32.5*	32.4 -32.7
	JJA	16	+38.5*	38.2 -38.8
	SON	9	+13.6*	13.5 -13.8

Table 2 .

 2 Same as in Table1but for actual evapotranspiration (AET).

	Period	Actual evapotranspiration (AET)
		% of stations	ΔAET (mm)	95 % HDI
	year	9	+45.7*	45.4 -46.1
	DJF	37	+8.29*	8.19 -8.40
	MAM	54	+23.2*	23.0 -23.3
	JJA	4	+7.69*	7.43 -7.96
	SON	4	+5.35*	5.25 -5.46

Table 3 .

 3 Same as in Table1but for liquid precipitation (PRELIQ).

	Period	Liquid precipitation (PRELIQ)
		% of stations	ΔPRELIQ (mm)	95 % HDI
	year	0	+23.2*	22.0 -24.3
	DJF	1	-3.07*	-3.60 --2.56
	MAM	3	+0.28	-0.12 -0.69
	JJA	1	+3.53*	3.08 -3.99
	SON	0	+10.4*	9.90 -10.9