Maria Brodovoi

Kilian Gruel

Lucas Chapuis

Aurélien Masseboeuf

Cécile Marcelot

Martin Hÿtch

Frédéric Lorut

Christophe Gatel

In-Situ Electrical Biasing of Electrically Connected TEM Lamellae with Embedded Nanodevices

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In recent years, there has been an increased interest in understanding local behavior and physical phenomena of electronic devices at the nm-scale. One of the readily available and appropriate tools for studying a wide variety of materials and nanosystems is Transmission Electron Microscopy (TEM). Its broad sensitivity, coupled to progress in aberrationcorrection, allows modern TEMs to reach spatial resolutions less than 50 pm. Besides conventional morphology and local structural and chemical studies, TEM investigations based on off-axis electron holography can also determine the electromagnetic [START_REF] Gatel | Inhomogeneous spatial distribution of the magnetic transition in an iron-rhodium thin film[END_REF][START_REF] Andersen | Field tunable three-dimensional magnetic nanotextures in cobalt-nickel nanowires[END_REF][START_REF] Midgley | Electron tomography and holography in materials science[END_REF] and deformation fields [START_REF] Hÿtch | Nanoscale holographic interferometry for strain measurements in electronic devices[END_REF] at the nanometer scale. Even though analytical TEM routinely contributes to failure analysis of semiconductor devices [START_REF] Engelmann | Application of analytical TEM for failure analysis of semiconductor device structures[END_REF], improvement or development of new devices require a deeper understanding of local properties when in operation. Operando TEM experiments, i.e., in-situ TEM biasing of electronic devices, are however rarely attempted [START_REF] Twitchett | Quantitative electron holography of biased semiconductor devices[END_REF][START_REF] Ikarashi | In-situ electron holography of surface potential response to gate voltage application in a sub-30-nm gate-length metaloxide-semiconductor field-effect transistor[END_REF][START_REF] De Knoop | Determining the work function of a carbon-cone cold-field emitter by in situ electron holograph[END_REF][START_REF] Einsle | In situ electron holography of the dynamic magnetic field emanating from a hard-disk drive writer[END_REF]. On the other hand, having access to local electric fields and electric charge distribution mapped at the nanometer scale across a single device in operation would address critical information for future developments or technology design improvements and optimization. The resulting information can also be beneficial to failure analysis, where very precise high-resolution local observations of nanodevices in operation could be used to observe isolated defects or abnormalities in a single region. This paper proposes a methodology to study electronic devices under real operation conditions by coupling advanced focusedion beam (FIB) sample preparation and off-axis electron holography. Additional finite element modeling (FEM) will allow the experimental results to be compared to numerical simulations including factors like specimen geometry, FIBinduced damage, stray fields, and charging effect from electron beam interactions.

Experimental Sample preparation

Experiments were carried out on a MIM (Metal-Insulator-Metal) matrix of unitary nanocapacitors manufactured by STMicroelectronics via their 28 nm FD-SOI (Fully Depleted Silicon On Insulator) process. A capacitor is a passive electrical device used for energy storage and delivery in analog and RF circuit applications. It is embedded in the intermediate backend dielectric layers and is composed of a thin layer of dielectric material sandwiched between two electrically separated metal electrodes. Our device consisted of a common 35 nm titanium nitride (TiN) bottom electrode, a thin 11.5 nm tantalum pentoxide (Ta2O5) insulator and 55 nm TiN top electrode (Figure 1).

A specific sample preparation, combining a TEM lamella preparation flow and a FIB-Edit capability, is required to bring electrical stimuli to the device while ensuring its electrical functionality and compatibility with TEM observations. This step is challenging, as the approach and methods for circuit modification and electrical connections must be adapted to the sub-µ scale. Samples are lifted-out from a fully processed wafer and prepared for electron transparency and in-situ biasing using a ThermoFisher Scientific TM Helios NanoLab 1200AT DualBeam TM with a Gallium liquid-metal ion source. The equipment includes a multiple line gas injection delivery system (MultiChem) and a gas injection system (GIS), both allowing ion and electron beam-induced deposition and selective etching for effective and high-rate circuit modification.

A conductive surface is essential for optimal charge evacuation during sample navigation and preparation; therefore, a thin conductive layer (gold-palladium alloy) is deposited on the specimen's surface, before any further steps, to ensure charge evacuation. The sample is then grounded by opening the pad (Fig. 2.a) and connecting the top surface to the grounded Si substrate: the charge from scanning electron beam (SEM) will be drained off to the ground during sample preparation.

Capacitors are susceptible to electrostatic discharge (ESD) damage; therefore, it is important to prevent the structure from floating during FIB-Edit. This can be done by connecting, in advance, both top and bottom electrodes. This is achieved by milling a line on the top surface and filling it with Pt (Fig. 2.a). This protective layer serves multiple purposes: it will act like a shield and preserve the structure during FIB milling, it will prevent charge accumulation, and it will temporarily bring both electrodes to an equipotential. The sample is then top-down deprocessed using FIB until the layer of interest is reached. Cu lines which connect the nanocapacitors are now exposed (Fig. The extracted sample is directly welded to an electrontransparent Si3N4 grid with Au prepatterned tracks (Fig. 2.e,f). W depositions are used to ensure electrical connection between the sample and grid. To meet electron holography criteria, the region of interest (here, Ta2O5 oxide) must be no more than 500 nm from the vacuum region (Fig. 2.g) where the reference wave will be. The next preparation step consists of separating the top electrode from the bottom one to create two distinct electrical connections for in-situ biasing. Two separation cuts (bottom left and top right) are visible on Fig. 2

.h.

A final thinning of the lamella with a low energy ion beam (8 or 16 kV) for minimizing surface damage is performed until a uniform thickness below 100 nm is reached for electron transparency.

TEM experimental setup

To perform in-situ electrical biasing experiments, a dedicated experimental setup must be implemented that consists of several levels of electrical connections (Fig. 3). These are required to bring electrical stimuli from the external power supply to the device embedded in the TEM lamella. The setup has been optimized to avoid electrostatic discharge (ESD) damage when inserting and observing the prepared device. As described previously, the TEM lamella (Fig. 3.e) is already welded to a grid with electrical contacts (Fig. 3.d) to bias the nanocapacitors. The Si3N4 grid is then mounted on a single-tilt Gatan HC3500 TEM holder (Fig. 3.b) equipped with 4 electrical contacts to input/output electrical signals from an external current/voltage generator. The grid is electrically connected to the holder contacts with thin conductive wires via µ-bonding (Fig. 3.c). Before operando TEM experiments, electrical testing (resistance measurements) is performed on the final connected electron transparent lamella to confirm that the nanocapacitors have not been electrically shorted. It should be pointed out that

a) b) c) d) e) f) g) h)
a resistance of several hundreds of ohms is acceptable for subsequent electrical biasing. TEM experiments were carried out on the I2TEM-Toulouse microscope, an HF3300-C (Hitachi) TEM specially designed for in-situ electron interferometry experiments. The microscope is equipped with a cold field-emission gun (CFEG) for optimal brightness and coherence, a double stage configuration and 2 post-specimen biprisms. The latest generation of imaging aberration corrector (BCOR from CEOS) allows for on and offaxis aberrations to be corrected in either Lorentz or conventional stage mode allowing wide fields of view and spatial resolution of 0.5 nm for the Lorentz mode and 87 pm in conventional HRTEM imaging [START_REF] Snoeck | Off-axial aberration correction using a B-COR for Lorentz and HREM Modes[END_REF]. Completing the setup is a high-speed 4K by 4K digital camera (OneView, Gatan Inc.).

Electron holography principle

The principle of off-axis electron holography relies on interfering two parts of a highly coherent electron: the object wave, which interacts with the sample and the surrounding electromagnetic fields, and the reference wave that has not undergone interaction with any field (Figure 4). The resulting interference pattern, called hologram, encodes the phase shift that the electron wave experienced when interacting with the local electromagnetic potentials, mandatory for quantitative mapping of electric fields inside and outside the device [START_REF] Tonomura | Applications of electron holography[END_REF][START_REF] Mccartney | Quantitative measurement of nanoscale electrostatic potentials and charges using off-axis electron holography: Developments and opportunities[END_REF].

The phase change caused by the electrical biasing and recorded in the x-y plane of the hologram is related to the electric potentials in the following way: 𝜙𝜙(𝑥𝑥, 𝑦𝑦) = 𝑐𝑐 𝐸𝐸 ∫ 𝑉𝑉(𝑥𝑥, 𝑦𝑦, 𝑧𝑧)𝑑𝑑𝑧𝑧 (1) where cE is a constant depending on the accelerating voltage of the fast electron and z the direction of the electron path. V is the sum of two contributions:

𝑉𝑉 = 𝑉𝑉 𝑀𝑀𝑀𝑀𝑀𝑀 + 𝑉𝑉 𝐸𝐸 (2)
with 𝑉𝑉 𝑀𝑀𝑀𝑀𝑀𝑀 related to the mean inner potential (MIP) proportional to the atomic number of the encountered atoms along the electron path and 𝑉𝑉 𝐸𝐸 the electric potential due to the applied bias. As a consequence, the measured phase shift can be separated in two contributions 𝜙𝜙 𝑀𝑀𝑀𝑀𝑀𝑀 and 𝜙𝜙 𝐸𝐸 . Only the latter is of interest in our study and will have to be extracted from the measured phase shift. Fortunately, 𝜙𝜙 𝑀𝑀𝑀𝑀𝑀𝑀 is independent of the applied bias and can be obtained at 0V before being subtracted from the total phase. The integration along the path of the incident electron means that the phase is sensitive to the electric potential within and external to the specimen (i.e., includes the contribution from stray fields).

Results

Electron holography experiments

Operando electron holography experiments were performed with an operating voltage of 300 kV, Lorentz stage, elliptical illumination and 2 post-specimen biprisms for flexibility in the holographic configurations and to eliminate the Fresnel fringe artefacts [START_REF] Harada | Double-biprism electron interferometry[END_REF].

Holograms were acquired during a long exposure time of 120s using dynamic automation [START_REF] Gatel | Unlimited acquisition time in electron holography by automated feedback control of transmission electron microscope[END_REF] and π-shift method [START_REF] Volkov | Double resolution electron holography with simple Fourier transform of fringe-shifted holograms[END_REF]. The fringe spacing was set to 1.1 nm (7 pixels). On-line phase images were displayed with the software package, HoloLive! (HREM Research Inc.), a plug-in for Digital Micrograph (Gatan Inc.). Off-line data processing was carried out with a 5 th order Butterworth Fourier-space filter centered on the sideband using routines based on Fourier analysis similar to those developed for the processing of dark-field electron holograms HoloDark 2.0 (HREM Research Inc.) [START_REF]HoloDark for DigitalMicrograph[END_REF]. The numerical mask size was chosen to obtain a spatial resolution of 2 nm in the final phase images. The projector and camera distortion-induced phase modulations were removed using a reference hologram

a) b) c) d) e)
recorded in the field-free vacuum area. The nanocapacitor studied is shown Fig. 5. For each applied bias, 3 different holograms were recorded for reconstructing the whole field of view. Fig. 6.a and 6.b correspond respectively to the amplitude image and the reference phase image at 0V (𝜙𝜙 𝑀𝑀𝑀𝑀𝑀𝑀 ) extracted from the holograms. Fig. 6.c and 6.d show the induced electric phase shift at 0.5V and 1V of applied bias respectively after subtracting the reference phase image at 0V: they only contain the electrostatic contribution of excess charges (due to voltage between the electrodes). As expected, the phase shift appears uniform inside the electrodes while a stray field is visible at the edges of the upper electrode due to its finite lateral dimension. The local phase variation located in the Ta2O5 oxide layer is associated to the electric voltage between the two TiN metal electrodes and consequently to the presence of a strong electric field (Fig. 6.d).

Three-dimensional (3D) finite element modeling (FEM) was performed with multiphysics software from COMSOL Multiphysics® (v.5.2 and license N°1022241) to have a deeper understanding and quantification of the electric potential distribution as a function of applied bias. The simulations enabled us to assess the impact of the 3D sample geometry and artefacts of preparation, namely FIB-induced sample degradation, on the final device behavior. The FIB-induced preparation artefacts to be considered are the local modification of the conductivity (i.e., conductive damage layers on both sides of the TEM lamella) and the creation of point defects which may modify the electric state of the sample. In addition, modeling allows us to separate the effects of the internal and external/leakage fields.

Regions of different mesh sizes were used to improve the efficiency of the calculations to create a model that sufficiently encompassed the device and surrounding vacuum. The model was surrounded by a layer that approximates an extension of the equations to infinity. Then, the mesh evolves into a free triangular type between the external boundaries and the nanocapacitor. The maximum mesh size in this area is 0.4 µm. The mesh dynamically adapts in size to meet the effective area of the nanocapacitor with a finer mesh size of about 2 nm minimum. The model is further divided into two main parts representing two distinct physics and allowing the simulation of conductive and insulating regions: Ohm's law (1) for induction and Gauss' law (2) for electrostatics, respectively. To achieve this, the variable V solved by these two parts must be 𝜎𝜎 : conductivity S/m ; 𝐽𝐽 𝑒𝑒 �� �⃗ : external (injected) current density A/m 2 ; 𝑄𝑄 𝑗𝑗 : current source A/m 3 .

-∇ • �𝜀𝜀 0 ∇𝑉𝑉 -𝑃𝑃 �⃗ � = 𝜌𝜌 (2) 𝜀𝜀 0 : vacuum permittivity F/m ; 𝑃𝑃 �⃗ : electric polarization vector C/m 2 ; 𝜌𝜌 : space charge density C/m 3 . Fig. 7 shows the comparison between the experimental (Fig. 7.b) result obtained for an applied bias of 1V and simulated (Fig. 7.d) phase image. We can see the excellent quantitative agreement, confirmed by the phase profiles extracted through the electrodes and the insulating Ta2O5 layer (Fig. 7.c). Simulations confirm the important linear variation of the phase shift which occurs in the oxide layer whilst the curvature of the phase in the electrodes caused by the stray field around the sample (see Fig. 4) is well reproduced. However, this agreement has been obtained for a reduced potential drop between the electrodes.

Since the sample was thinned with a Ga-based focused ion beam, the FEM model also includes conductive layers around the sample. This reduces the potential drop across the dielectric with respect to the applied voltage due to current leakage in the conductive layers. Parameters such as the thickness (due to amorphization on the edges of specimen) and the conductivity of the surface layers were adjusted to produce the best fit with the experimentally measured phase shift (Fig. 7.c). Experimental and theoretical curves may however present discrepancies (like, notable curvature inside the dielectric, related to the evolution of accumulated charge in the oxide) which may need further adjustments by introducing different conditions to the finite element model.

Conclusions

A set of MIM nanocapacitors have been extracted from production lines and biased for in-situ electrical experiments within the TEM. With the device functioning in real working conditions, operando electron holography allowed the electric potential to be mapped across the active areas and the surrounding stray fields. By comparing the experimental results with FEM modeling, local electrical parameters could be deduced such as electrostatic potential and electric field inside the dielectric layer, capacitance, and surface charge density. The study represents a complete in-situ electrical investigation of a single device under fixed bias and further interpretation of these results will help to enhance sample preparation and minimize its impact on physical parameters modification.

There is no reason that the approach cannot be applied to more complex electronic devices (transistor, phase change memory) by modifying the sample preparation workflow to achieve a successful circuit modification and enable electrical stimulation of the device. We hope that the ability to study local electric fields in targeted devices from production lines will help efforts to improve devices in terms of reliability, speed, and power consumption.

Figure 1 :

 1 Figure 1: TEM image of MIM nanocapacitor stack.

  2.b) and once again, top and bottom electrodes are brought to the same potential by depositing a Pt layer (Fig. 2.c).

Figure 2 :

 2 Figure 2: a) Sample grounding & electrode connections. b) Visible Cu lines after FIB deprocessing. c) Protective Pt layer. d) Sample before lift-out. e) Sample welding to the dedicated TEM grid. f) Chunk welded on the TEM grid. g) TEM grid membrane milling & free area for EH reference wave. h) Final thinning & electrodes separation.

Figure 3 :

 3 Figure 3: Experimental setup; electrical connectivity levels: ab) TEM microscope -sample holder for electrical biasing via external electrical wiring. b-c) Sample holder -specific TEM grid via µ-bonding. d-e) TEM grid -device embedded in the TEM lamella via FIB-assisted conductive layers.

Figure 4 :

 4 Figure 4: Principle of electron holography and link between potential distribution made by a capacitor and the observed phase.

Figure 5 :

 5 Figure 5: Low magnification TEM micrograph of the studied MIM nanocapacitor. The red dotted line represents the position of the biprism.

Figure 6 :

 6 Figure 6: a) Amplitude image of the capacitor C. b) Phase image under 0V bias used to remove the mean inner potential contribution. c) Phase image of the induced electrostatic potential under 0.5V bias. d) Phase image of the induced electrostatic potential under 1V bias.

Figure 7 :

 7 Figure 7: Comparison of experimental and simulated data a) Amplitude image of the studied capacitor. b) Phase image of the same capacitor biased under 1V. c) Profiles of the experimental and simulated phase image extracted along the black arrows on the phase images. d) Simulated phase image from FEM.

Acknowledgments

The authors acknowledge the French National Research Agency under the "Investissement d'Avenir" program reference No. ANR-10-EQPX-38-01" and No. 11-IDEX-0002, the "Conseil Regional Midi-Pyrénées" and the European FEDER for financial support within the CPER program. This work was also supported by the French national project IODA (ANR-17-CE24-0047). This project (ADVENT -16ENG06) has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme. The research leading to these results has received funding from the European Union Horizon 2020 research and innovation programme under grant agreement No. 823717 -ESTEEM3.


Yong He
email: yong.he@uca.fr

An Empirical Test on Regional Spillovers through Intra-and Inter-national Trade

Keywords: intra-national trade, intra-national spillovers, regional input-output tables, regional disparity. JEL Classification: F1, R1, O1, O3, O4

This study aims at providing econometrical evidence to support a theory predicting that in developing countries, intra-national trade, together with international trade, form a network in which high growth achieved by the developed regions spills over to the less developed regions. Using China's 2007 foreign trade data and provincial input-output tables, the key variables on intra-and international imports of technological inputs are built for estimating their impacts over outputs at the province and sector levels. It is found that in the less developed regions, intra-national imports rather than international ones made significant contribution to production. In the developed regions, these impacts were just inversed. This result confirms the theory on the existence of the trade network.

Introduction

Within a developing country with huge regional inequality, intra-national trade must play a decisive role for regional development. It can be observed that the developed regions have stronger trade connections with foreign countries, while the less In this study, we first derive an econometrically testable equation of production function including distinct intra-national and international imports of technological inputs. Then, we construct the variable of international imports of equipment inputs by province and sector based on foreign trade data, and the variable of intra-national imports of equipment input by province and sector based on 2007 provincial inputoutput tables. As expected, the estimation results strongly support the theoretical predictions.

This paper is organized as follows: following this introduction, section 2 introduces the background from which the theoretical model is derived and the main predictions of the model to be tested. Section 3 exposes the estimation method. In Section 4, I present the data and the construction of the variables. Section 5 analyzes the results before concluding.

Theoretical Background and Issue

This section introduces the main intellectual sources that inspire and support the construction of the theoretical framework on spillovers of intra-national trade, which constitutes the basis for our empirical test. We first present the two types of intranational trade: the one driven by comparative advantage, and the other driven by increasing returns and agglomeration among regions, and then provide arguments on why the former, albeit its limitations, is more appropriate for a developing economy. Afterward, we explain the necessity for choosing a model of technological spillovers with an expanding variety of intermediates to describe regional development through intra-national trade in developing countries. In the end, we introduce the issue for testing.

As [START_REF] Krugman | Interregional and International Trade: Different Causes, Different Trends? in Regional Science Matters[END_REF] puts it, there are two driving forces to intra-national trade: the comparative advantage, and increasing returns and agglomeration. On the one side, He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 4 © Jadavpur University [START_REF] Ohlin | Interregional and International Trade[END_REF] on interregional trade is a more general approach that covers at once intra-and inter-national trade. In Ohlin's view, "regions have different factor endowments, while the factors within a region are essentially similar. …it is assumed that the factors of production are inter-regionally immobile but intra-regionally freely mobile." (op. cit., p.5) Each region specializes in the productions for which its factors are relatively abundant and cheap, and then all regions gain from trade than staying in autarchy.

One the other side, a region could be defined as a district à la Tünen in which individuals and firms relocate to where increasing returns are higher [START_REF] Isard | Location and Space-economy[END_REF]). In new economic geography, both interregional and international trade are described from the point of view of location theory and agglomeration.

The explanatory power of the new economic geography depends crucially on the extent of possible factor mobility and an agglomeration process to take place. In a developed economy, as factor mobility is high and agglomeration process occurs well, regional disparity becomes endogenous. Exogenous regional disparity refers to what [START_REF] Cronon | Nature's Metropolis, Chicago and the Great West[END_REF] distinguishes endogenous inequality from the exogenous one, which he also calls first nature inequalities; there are natural differences in resource endowment, climate feature, and geographical characteristics among regions. In contrast, endogenous, or second-nature inequalities, are the result of human actions to improve upon first-nature inequalities. Combes et al. (2008, p15) claim that the new economic geography is more appropriate to approach to endogenous regional disparity for the developed world.

When applied to developing countries, this approach exposes serious limitations, because in these countries, the industrialization process is relatively new. Conglomeration, which parallels this process, is not yet a major force. Due to the high costs of mobility associated with the backwardness in transport infrastructure and institutional impediments, factor movement, labor, and capital meet serious obstacles.

Ohlin's comparative advantage approach is more apt to characterize the intranational trade in developing countries, because one of the intrinsic features that distinguish a developing country from a developed one is that while in the former, regional disparity remains wide, in the latter, agglomeration and factor relocation He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 5 © Jadavpur University have made it much smaller. Using Krugman's (2015) expression, "America is flat" and "Americans are doing pretty much the same thing everywhere". In developing countries, however, in the absence of consequential agglomeration forces, exogenous regional disparity keeps dominant. This distinction is supported by some well-known studies such as [START_REF] Williamson | Regional Inequality and the Process of National Development: a Description of the Patterns[END_REF] which shows that during the development stages of a country, regional disparity tends to be larger and will increase up to a certain point. After this point, it will decrease due to labor mobility that reduces income inequality between regions. More recently, Lessmann (2014) makes use of a unique panel data set of spatial inequalities in 55 countries at different stages of economic development, covering the period 1980-2009, and provides strong support for the existence of an inverted U relationship of within-country regional inequality. Whereas Ohlin's view accommodates better developing countries, it has a major limitation in dealing with the issue of technological spillovers through intra-national trade. Spillovers on the basis of regional differences in factor endowment are difficultly modeled. On the contrary, as in most developing countries, regional disparity is also in accordance with differences in technology, always keeping the spirit of Ohlin's comparative advantage-based intra-national trade, there is a need to assign it a new form: a technological spillover model.

A vast literature on international technological spillovers provides a spur for our approach [START_REF] Findlay | Relative Less Developedness, Direct Foreign Investment, and the Transfer of Technology: A Simple Dynamic Model[END_REF]; [START_REF] Krugman | A Model of Innovation, Technology Transfer, and the World Distribution of Income[END_REF]; [START_REF] Dollar | Technological Innovation, Capital Mobility, and the product Cycle in North-south Trade[END_REF]; [START_REF] Grossman | Trade, Knowledge Spillovers and Growth[END_REF]). To extend these models to regional analysis, the marked regional disparity within developing countries implies that, first, there is a strong internal technological demand and, second, the less developed regions have a relatively low technological absorption capability. The most important argument justifying intra-national spillovers is technological distance. [START_REF] Fu | International and Intra-national Technology Spillovers and Technology Development Paths in Developing Countries: the Case of China[END_REF] have shown that within developing countries, regional technology transfer takes place more effectively when technological distance is small. In connection with this argument, relating to international trade, intra-national trade can be conjectured to offer lower but more appropriate technology and to be less costly in transport and more efficient in technology transfer. This is due to the proximity in language, institutions and culture, and facilities in labor and capital mobility, knowledge communication, and learning and formation. He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 6 © Jadavpur University Two kinds of models have been generally applied to delineate technological spillovers: those with an expanding variety of intermediates, and those of quality ladders. We choose the first approach (grounded on the seminal work by [START_REF] Ethier | National and International Returns to Scale in the Modern Theory of International Trade[END_REF] and [START_REF] Romer | Endogenous Technological Change[END_REF]) for its facility of adaptation to a trade model. Owing to the use of higher technological input, the final products are upgraded in quality, in design, and in variety even without necessarily buying sophistical equipment or changing production process. A stream of econometric papers such as Feenstra et Grounded on all these theoretical developments, in He (2017), a dynamical model is built to analyze the nature of intra-national trade in developing countries. There is a developing country, within which there are multiple regions ranked by technological capability and characterized by their differences in number of varieties of intermediate goods at the starting time. There is also a developed foreign country with a technological level higher than all regions in the developing country. A region,as its production is a function of the number of varieties of intermediates, must choose between inventing, or importing and then imitating new intermediate goods. If importing, a region must also choose how much to import, and then importing from where: other advanced regions or a developed country. Intra-national trade shares of the regions are shown to be endogenously determined in the light of a theory of choice between intra-and international trades. This choice is contingent on the differentials in trade costs and in assimilation costs between intra-and international imports, with the focus on high assimilation costs of new intermediates from foreign trade for the less developed regions. In our approach, the one-to-one dynamical analysis is extended to multiple regions in which their choices on the basis of different technological capabilities lead to the formation of a network connecting inter-and intra-national trades. Within the network, intra-national trade network prolongs the technological spillover effects emanating from international trade, through the advanced regions, to the least developed regions, so that all the partners involved in the trade network, originally with various growth rates, reach the highest growth rate benefited initially by the foreign country.

The main result of the model is that regional spillovers occur within the trade network in which regions choose to import technological inputs from the regions (or He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 7 © Jadavpur University the foreign country) with which they have narrower technological distance. The number of layers depends on the degree of regional inequality, geographical, and other factors. At any layer of the trade network, the optimal number of varieties of intermediates will be altered by that of the highest layer. Even though in the beginning, the less developed regions have low growth rates, the highest growth rate benefitted by the foreign country will be transmitted to all levels of the trade network.

The above theoretical results give rise to two empirically testable predictions: Prediction 1: In a developing country, the developed regions are the main international importers of technological inputs, and intra-national import ratios of technological inputs are much higher in the less developed regions than in the developed regions.

Prediction 2: The developed regions mainly benefit from spillovers through international imports of technological inputs. Inversely, the less developed regions mainly benefit from spillovers coming from intra-national imports of technological inputs.

For confirming the Prediction 1, in He (2017, section IV.B), several tables with descriptive statistics on the differences in share of intra-national trade over total trade were made. The main drawback of this study is the lack of an econometrical test on the regional technological spillover effects to validate Prediction 2.

The empirical tests of spillovers via intra-as well as international imports of technological inputs remain a challenge because they meet at once methodological obstacles and data availability. Testing intra-national spillover effects, in a strict sense, requires panel data with multiple periods and various trade levels. The best way to test the spillover effects will be employing the panel data covering a period enough long and a number of regions of different development levels. Nevertheless, intra-national trade data at the regional level are scarce. The time series data on intranational trade are much scarcer. In the face of this constraint, we choose an alternative way that does not require time series database.

Given the objective is to find the evidence of the spillovers via intra-national trade for the less developed regions, and the spillovers via international trade for the more He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 8 © Jadavpur University advanced regions, our econometric test is guided by the following reasoning. If among the less developed regions, a higher share of intra-national imports of equipment allows some regions improving their production relating to the others keeping a low share of these imports, this is a proof of the existence of intra-national spillovers. The same reasoning applies to the international imports of equipment by the advanced regions. This result could confirm the Prediction 2.

Therefore, in the following sections, we first derive an econometrically testable equation of production function including distinct intra-national and international imports of technological inputs. Then, for our test, we construct the variable of international imports of equipment inputs by province and sector based on foreign trade data, and the variable of intra-national imports of equipment input by province and sector based on 2007 provincial input-output tables.

Estimation Method

As a large country and having experienced lasting economic growth and technological progress, China is an appropriate case for evaluating the effects of intra-national trade for regional development. China's 31 provinces are conventionally classified into three large regions: Costal region (10 provinces), Central region (9 provinces), and Western region (12 provinces, with Tibet ruled out in this study due to data missing). Coastal region is the most developed and the main exporter of China to the world. It is followed by Central region, and Western region is the less developed. Their population shares were 36.6%, 35.6% and 27.8%, and their GDP shares 55.3%, 27.4% and 17.3% in 2007. The key idea guiding our empirical tests is to check the differences in effects of international and intra-national imports of technological inputs over output in three regions, respectively. In other words, total sample will be divided into three sub-samples, and the regressions will be made with them separately.

The production function of the province h's output of sector i takes the following form: Ea refer to international imports and intra- national imports of technological input, respectively.  is a parameter reflecting the extent to which international imports of technological input is different to intranational imports of technological input in technological content. If 1   , they make equal contribution to output. But in accordance with our theoretical framework,  for the developed regions could be larger than for the less developed regions.

x l k ih ih ih ih ih Ea Ee L AK Y     ] [   (1)
In logarithm form, Equation ( 1) becomes

] [ ih ih x ih i ih k ih Ma Me ln lnL lnK A ln Y ln          (2)
The last term is dealt with using linear expansion around

i ih a Me  and i ih b Ma  in which i
a is set as the mean value of international imports of technological input by sector across 30 provinces, and i b is the mean value of intra-national imports of technological input by sector across 30 provinces. This term becomes

i i ih i i i x i i ih i i i x i i x b a Ma b a b b a Me b a b b ln(a                ] 1 ) [ (3)
Collecting the first term into the error term, we get the equation for testing: Intuitively, to be able to affirm positive trade spillover effects, two outcomes are required. 1) Intra-national imports of technological input have a statistically significant positive impact to the output for the less developed regions whereas they are insignificant for the developed regions; 2) International imports of technological He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 10 © Jadavpur University input have a statistically significant impact to the output for the developed regions whereas they are insignificant for the less developed regions. These results could confirm Prediction 2.

ih i i ih xa i i ih xe ih i ih k ih b a Ma b a Me lnL lnK A ln Y ln              (4)
If both variables are positively significant for both regions, the comparison of the levels of the coefficients will be needed. To be in accordance with Prediction 2, the coefficient of international imports of technological input must be larger for the developed regions than for the less developed regions, and the coefficient of intranational imports of technological input must be larger for the less developed regions than for the developed regions.

Data and Construction of the Variables

In total, there are input-output tables of 30 provinces with 42 sectors: 26 industrial sectors (including 17 manufacturing sectors, in which there are 6 equipment sectors), 15service sectors, plus agricultural sector. 1Based on these provincial 2007 input-output tables and 2008 statistic yearbooks (containing 2007 data), we provide some descriptive statistics on Chinas inter-and intra-national trade and their distribution among the three grand regions. In 2007, China's intra-national trade (defined as inter-provincial trade volume, without taking into account of intra-provincial trade) account to 72.69% of the GDP, more than 2.21 times of foreign trade. Intra-national imports of manufactured products account to 1.91 times of international imports of manufactured products. Intra-national imports of equipment products were 1.42 times of international imports of equipment products. These numbers indicate the crucial importance of intra-national trade for such a country like China, albeit with high international openness.

The apportionment of inter-and intra-national imports is highly unequal. On international manufacturing imports, the shares were 91.09%, 6.13%, and 2.78% for the Coastal, Central and Western regions, respectively. Their shares on intra-national He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 11 © Jadavpur University manufacturing imports were 55.34%, 26.38% and 18.28%, respectively. These numbers indicated that while foreign trade is largely concentrated in coastal region (due to the central role of Chinese coastal region as "world factory"), intra-national trade is relatively more important for the less developed regions. The less developed regions have larger dependence on intra-national equipment imports. The net intranational imports of equipment were -51.7, 91.4 and 87.6 billion USD for the three regions, respectively. Only the Coastal region is net exporter to other regions while the less developed regions are net importers.

Based on these provincial input-output tables and statistic yearbooks, we built all the variables for testing Equation (4) defined in Table 1. Afterward we introduce one by one their construction.

Table 1 The definitions of the variables

Variable name Definition

Lnoutput

The output by province and sector in logarithm form.

Lncapital

Returns to capital by province and sector in logarithm form are used as the proxy of capital.

Lnlabor

Compensation of employees by province and sector in logarithm form is used as proxy of labor.

Rate_inter_equip

The share of international imports of equipment by province and sector in i i b a  defined in note (4) of this table.

Rate_intra_equip

The share of intra-national imports of equipment by province and sector in Notes: 1) international imports of equipment goods are calculated based on 2008statistic yearbooks of the provinces; 2) intra-national imports of equipment goods are computed on the basis of "inflow" and "outflow" by province and sector in China's 2007 provincial input-output tables and of international trade data in the 2008 statistic yearbooks of the provinces; 3) intra-national imports are specified as inter-provincial imports without intra-provincial imports included; 4) is set as the mean value of international imports of equipment input by sector across 30 provinces, and is the mean value of intra-national imports of equipment input by sector across 30 provinces.

Lnoutput is based on the output by province and sector that are found in all inputoutput tables. To build Lncapital, capital is surrogated by capital returns that are composed in three terms in input-output tables: net taxes on production, operating He / Regional Spillovers through Trade

Trade and Development Review, Vol. 12, Issue 1-2, 2019 12 © Jadavpur University surplus, depreciation of fixed capital. For Lnlabor, labor income is reflected by the item: compensation of employees is used as labor inputs. Usually the amount of labor hired and the book value of capital are used as labor and capital variables. These data are, however, absent in input-output tables. Here assuming that labor and capital inputs are valued according to their marginal productivities (capital price is the present values of the net cash flows), these items provide a convincing measurement of these inputs. One advantage of this method is that these values incorporating the differences in wage and capital returns as proxies reflect better regional disparity on capital and labor, because, labor amount in quantitative terms are comparable, but highly unequal in qualitative, or human capital term among the regions.

In what follows, the issue is how to construct the two last explanatory variables. The first Chinese provincial input-output tables were published in 2002. In 2007, for the first time, in them, there are "inflow" (international imports plus inter-provincial imports) and "outflow" (international exports plus inter-provincial exports) by province and sector. With the data on international trade 2007 by sector and province, which can be found in the statistic yearbooks 2008 of these provinces, "inflow" is potentially usable for constructing the inter-and intra-national imports variables. 2The first possibility was using the data on manufactured inputs, and testing the impacts of inter-and intra-national imports of manufactured goods on regions' outputs. This method, however, meets a serious concern on the bias coming from the fact that China's manufactured imports are to very large extent for the purpose of processing. 3 As the consequence of the processing, a large volume of inter-and intra-national imports of manufactured goods are driven by the re-exportation of their final goods. This is not appropriate for testing the impact of the imports He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 13 © Jadavpur University motivated by enhancing technological capabilities. Instead, we choose to make our test with the data of equipment goods. We reason that as a developing country, China's imports of equipment goods have a limited processing nature, and a substantial share of them is for fitting their technological gaps.

Based on available data on equipment imports from the statistic yearbooks 2008 of 30 provinces and of the equipment inflows from the input-output tables 2007 of 30 provinces, we must find some methods to apportion the intra-and internationally imported equipment goods among sectors within the same province. We leave the detailed presentation on methods of estimation in the Appendix. With these estimations, we get the Rate_inter_equip and Rate_intra_equip defined in Equations ( 3) and ( 4).

With 30 provinces of 42 sectors, in total 1237 observations are obtained (the sector "Scrap and waste" and exceptionally some other sectors had missing values, and 23 were dropped). Instead of testing the impact on 42 sectors, we merely test the impacts on 26 industrial sectors and 17 manufacturing sectors. In other 15 service sectors and the agricultural sector, the imported equipment inputs are either absent or very weak. Therefore, it would not be pertinent to test the impact of imported equipment in these sectors.

In Table 2, only the descriptive statistics of industrial sectors are presented. It can be observed through Rate_inter_equip and Rate_intra_equip that albeit in absolute terms, Costal region employed at once more inter-and intra-national imports of equipment goods, in relative terms, Coastal region worked with more international imports than intra-national imports of equipment inputs. Western region employed much less international imports than intra-national imports of equipment inputs. Central region kept the intermediate level of them.

Our one-year data is a panel of two dimensions: province and sector. Applying pooled ordinary least squares to panel data might be overly restrictive and can have a complicated error process (e.g., heteroskedasticity across panel units, and (or) serial correlation within panel units). For this reason, panel-data estimation method is employed. Sectoral fixed effects are important, while as provinces are clustered into three grand regions according to their development levels, within one region, province fixed effects must be less important. Therefore, group variable is set as He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 14 © Jadavpur University sector. On the choice between the fixed-effects (FE) and random-effects (RE) models, Hausman tests favor the former. Therefore, the fixed-effects (within) regressions are chosen. Notes: 1) international imports of equipment goods are calculated based on 2008 statistic yearbooks of the provinces; 2) intra-national imports of equipment goods are computed on the basis of "inflow" and "outflow" by province and sector in China's 2007 provincial input-output tables and of international trade data in the 2008 statistic yearbooks of the provinces; 3) intra-national imports are specified as inter-provincial imports without intra-provincial imports included;4) the provinces are clustered in three regions according to conventional method with 10 provinces in Coastal region, 9 provinces in Central region and 11 provinces in Western region (with Tibet being ruled out); 5) the definitions of the variables are made in this section. Notes:1) international imports of equipment goods are calculated based on 2008 statistic yearbooks of the provinces; 2) intra-national imports of equipment goods are computed on the basis of "inflow" and "outflow" by province and sector in China's 2007 provincial input-output tables and of international trade data in the 2008 statistic yearbooks of the provinces; 3) intra-national imports are specified as inter-provincial imports without intra-provincial imports included;4) the provinces are clustered in three regions according to conventional method with 10 provinces in Coastal region, 9 provinces in Central region and 11 provinces in Western region (with Tibet being ruled out); 5) the definitions of the variables are made in section IV; First, the variables reflecting both intra-and international imports of equipment inputs have their signs and significances in accordance with the theoretical prediction: For Coastal region, international imports of equipment inputs exerted He / Regional Spillovers through Trade

Results and Analysis

Trade and Development Review, Vol. 12, Issue 1-2, 2019 16 © Jadavpur University significant positive impacts on production, whereas intra-national imports of equipment inputs were insignificant. 4 For Western region, intra-national imports of equipment inputs had significant positive effects. In contrast, international imports of equipment inputs were not positively significant. The impact of intra-national imports of equipment input was stronger in manufacturing sectors than in industrial sectors. For Central region, while based on 26 industrial sectors,the international imports of equipment inputs are significant at the5% level, based on 17 manufacturing sectors, this significance disappears, indicating that the impacts of these inputs were weaker for Central region than for Costal region. This is also in concordance with the theoretical prediction.

To summarize, the econometric results provide clear evidence that intra-national trade was more beneficial than international trade to the less developed regions in terms of the productive contribution whereas international trade was more beneficial than intra-national trade to the developed regions. Thus, with the validation of Prediction 2, we confirm the existence of a trade network within which there are technological spillovers from the developed to the less developed regions through both inter-and intra-national trade.

One concern on the validity of the above estimations is: should the presence of endogeneity be suspected? [START_REF] Olley | The Dynamics of Productivity in the Telecommunications Equipment Industry[END_REF] and [START_REF] Levinsohn | Estimating Production Functions Using Inputs to Control for Unobservables[END_REF] have extensively discussed the presence of simultaneity and endogeneity in the case of the measurement of the impacts of intermediates on productivity. If inputs are chosen based on productivity shocks, a province with a higher productivity shock may use more imported inputs. Another possible source of endogeneity is that international exports shock as unobservable variable in error term may be correlated with the interprovincial imports of intermediates. In both cases, one of the conditions for unbiased and consistent estimation is violated. To deal with the endogeneity problem, in most previous work on the measurements of the impacts of intermediate inputs on productivity, panel data with multiple years were used. Tow-period data were needed for testing Granger causality [START_REF] Kim | The Effect of Imports and Exports on Total Factor Productivity in Korea[END_REF]). More often GMM estimator and Proxy Estimator following [START_REF] Olley | The Dynamics of Productivity in the Telecommunications Equipment Industry[END_REF] and Levinsohn and © Jadavpur University [START_REF] Levinsohn | Estimating Production Functions Using Inputs to Control for Unobservables[END_REF] were employed to compare with OLS estimator. Here as there is only one-year data, are the estimations still valid? Three arguments can be offered in favor of our approach.

First, in 2007, economic growth rates among provinces in China are unusually synchronized, with these rates for Coastal, Central and Western regions were 14.7%, 14.4%, and 13.7% respectively. Even though the growth rates by sector among provinces were likely to be more variable, these variances, shaped by the GDP variances, might be quite moderated. Thus it can be thought that productivity shocks on interprovincial equipment imports, even existing, were weak. Therefore, the estimations with the province-level data, instead of plant data, make sense.

Second, about another source of endogeneity: international exportation is a variable that affects at once the output and the manufactured imports, in such main Chinese processing and exporting provinces as Guangdong and Shanghai, this endogeneity caused by international exports must be strong. Therefore, there is a concern if measuring with manufactured inputs. Measuring with equipment inputs as we did, however, suffers from limited processing effects and reflects to large extent the demand for technological enhancement.

Last, most studies that measure the output impacts of imported intermediates based on plant level data, with different estimators, cannot lead to conclude that the results with OLS estimator without tackling endogeneity were systematically under or overbiased. For instance, [START_REF] Halpern | Imports and Productivity. Federal Reserve Bank of New York[END_REF], employing all Hungarian manufacturing firms during 1992-2003, got productivity impact of imports of 16.9 percent with OLS estimator, and 17.7 percent with OP estimator following [START_REF] Olley | The Dynamics of Productivity in the Telecommunications Equipment Industry[END_REF]. [START_REF] Kasahara | Does the Use of Imported Intermediates Increase Productivity? Plant-level Evidence[END_REF] on the basis of 3598 Chilean manufacturing plants from 1979 to 1996, got productivity impact of imports of 9.6 percent with OLS, 5.8 percent with GMM system, and 14.33 percent with Proxy Estimator.

With the above arguments, it seems reasonable to conclude that the endogeneity is not a serious concern and the suspect that results are significantly biased could be ruled out in this study. He / Regional Spillovers through Trade Trade and Development Review, Vol. 12, Issue 1-2, 2019 18 © Jadavpur University

Conclusion

This study used China's 2007 input-output tables and foreign trade data at the province level to test the existence of spillover effects via intra-and international trade analyzed in a theoretical model. It is found that in the less developed regions, intra-national imports rather than internationals imports of equipment inputs made significant contribution to the outputs in the production function. In the developed regions, these impacts are just inversed. These results confirm the existence of a trade network in which mainly the less developed regions benefited from spillovers via intra-national trade, while the developed regions gained this benefit via international trade.

  of international and intra-national imports of technological inputs by province and sector in i i b a  , respectively.

  in note (4) of this table.

Table 2 Descriptive statistics of the variables (industrial sectors)

 2

			Coastal region	Central region	Western region
			Obs. 250		Obs. 228		Obs. 280
			mean	S.D	mean	S.D	mean	S.D
	Lnoutput (	)	15.291	2.011	14.736	1.535	13.471	1.919
	Lncapital (	)	13.566	1.973	13.150	1.578	11.980	1.915
	Lnlabor (	)	12.583	2.023	12.348	1.568	11.183	1.911
	Rate_inter_equip	0.638	0.534	0.271	0.272	0.082	0.183
	Rate_intra_equip	0.936	1.594	0.690	0.960	0.422	0.610

Table 3

 3 presents the regression results for three grand regions. The high values of Rsquared, and of Wald chi2 on the significance of the regression relationship validate the chosen regression model. The values of Rho (fraction of variance due to individual effect) are in general large, signifying that the individual effects of sectors are strong, and panel estimators are better than pooled estimators. With 26 industrial sectors, the observation numbers are 250, 228, 280 for three regions, respectively. With 17 manufacturing sectors, the observation numbers are 169, 151, 184 for three regions, respectively.

	He / Regional Spillovers through Trade

Table 3 Regression results

 3

			Industrial sectors			Manufacturing sectors
		Coastal region	Center region	Western region	Coastal region	Center region	Western region
		lnoutput	lnoutput	lnoutput	lnoutput	lnoutput	lnoutput
	lncapital	0.670***	0.550***	0.554***	0.687***	0.555***	0.555***
		(0.0359)	(0.0307)	(0.0931)	(0.0463)	(0.0376)	(0.125)
	lnlabor	0.286***	0.397***	0.340***	0.272***	0.398***	0.321**
		(0.0262)	(0.0339)	(0.104)	(0.0320)	(0.0429)	(0.136)
	Rate_inter_equip	0.174***	0.122**	0.0806	0.169***	0.0552	0.276
		(0.0330)	(0.0589)	(0.176)	(0.0336)	(0.0712)	(0.324)
	Rate_intra_equip	0.00686	0.00838	0.175***	0.00680	-0.0191	0.216***
		(0.0109)	(0.0274)	(0.0376)	(0.0114)	(0.0648)	(0.0661)
	_cons	2.482***	2.562***	2.960***	2.544***	2.618***	3.244***
		(0.252)	(0.307)	(0.540)	(0.312)	(0.356)	(0.760)
	Wald chi2(7)	2470.92	470.94	419.78	2361.59	459.05	178.48
	(prob<chi2)	(0.000)	(0.000)	(0.000)	(0.000)	(0.000)	(0.000)
	R-sq
	(within	0.9724	0.9359	0.9002	0.9816	0.9371	0.8703
	between	0.9471	0.9529	0.9583	0.9354	0.9774	0.9852
	overall)	0.9611	0.9474	0.9296	0.9693	0.9565	0.9246
	Rho	0.5762	0.5315	0.4007	0.5320	0.3318	0.2574
	N	250	228	280	169	151	184

The equipment sector includes: 1) manufacture of metal products;

2) manufacture of machinery; 3) manufacture of transport equipment; 4) manufacture of electric machinery and instrument; 5) manufacture of electronic and communication equipment; and 6) manufacture of instruments, meters and other measuring equipment.

Note that we are conscious of some unavoidable limitations of the input-output data. First, in comparison with the data at the firm level, their aggregations at the province level may bring about substantial loss in microeconomic value. Second, the trade data: "inflow" and "outflow" are collected without distinguishing their origins and destinations. Furthermore, the aggregated nature of the data is often associated with omission, subjective arbitration and smoothness. We judge acceptable, though, to use them for testing our theoretical propositions only, rather than, for example, applying the model for forecasting purpose.

According to China Statistic Yearbook 2008, the volume of processing trade represented 58% of the Chinese foreign trade of manufactured goods in 2007.

The coefficients of Rate_inter_equip and Rate_inter_equip cannot be interpreted in the same way relative to the coefficients of the other variables. According to Equations (2) and (3), the former is not comparable to the latter.

We believe that this method of estimation with more subdivided equipment sectors could give rise to a more faithful result. This method, however, could not be applied to the estimations of sectoral distribution of internationally imported equipment inputs due to the incomplete detailed information on equipment sector.

© Jadavpur University

Appendix The method of apportioning international and intra-national imports across sectors

The first task is to approximate the international equipment imports by sector and province, because we have these data by province, but not by sector. From the statistic yearbooks 2008 of 30 provinces, their international imports of equipment goods are gathered. For estimating their distribution as inputs among sectors, a weighting method must be found. First, only a share of these goods was employed as inputs and the other as final consumption. Thus the total international equipment imports time the ratio: the equipment input/(equipment input + final consumption of equipment) result in the internationally imported equipment input by province. Second, the obtained internationally imported equipment inputs are apportioned among the 42 sectors according to the shares of their equipment inputs in the total equipment inputs of the province. 5 5 The weighting method has also been generally used in the estimation of multi-regional trade relationship by official statistic bureaus (cf. e.g., National Information Center 2005 p.20). Based on provincial input-output tables, we are able to compute the equipment inflow by province. Its repartition by sector must also be estimated. First, using the ratio: the equipment input/(equipment input + final consumption of equipment) at the province level, the equipment inflow used as inputs is estimated. Then the following formula is employed to distribute the equipment inputs among the 42 sectors. There are in total 6 equipment sectors. For one of the 42 sector, h, its equipment inflow is

where,

hj weight) input ( is the sector h's share in the j inputs of the province. 6 Finally, intra-national imports of equipment inputs by sector are obtained by subtracting international equipment imports from their equipment inflow used as inputs.)

Based on these estimated internationally and intra-nationally imported equipment inputs by sector and province, we are able to calculate the Rate_inter_equip and Rate_intra_equip defined in Equation ( 4).


Alexandre Pohl

Zunli Lu

Wanyi Lu

Richard G Stockey

Maya Elrick

Menghan Li

André Desrochers

Yanan Shen

Ruliang He

Seth Finnegan

Andy Ridgwell

Vertical decoupling in Late Ordovician anoxia due to reorganization of ocean circulation

Geochemical redox proxies indicate that seafloor anoxia occurred during the latest Ordovician glacial maximum, coincident with the second pulse of the Late Ordovician mass extinction.

However, expanded anoxia in a glacial climate strikingly contrasts with the warming-associated Mesozoic anoxic events, and raises questions both as to the causal mechanism of ocean deoxygenation and its relationship with extinction. Here we firstly report iodine to calcium ratio (I/Ca) data that documents increased upper-ocean oxygenation despite the concurrent expansion of seafloor anoxia. We then resolve these apparently conflicting observations as well as their relationship to global climate, by means of a series of Earth System model simulations. Applying available Late Ordovician (Hirnantian) sea-surface temperature estimates from oxygen isotope studies as constraints, alongside our I/Ca data, leads us to identify a scenario in which Hirnantian glacial conditions permit both the spread of seafloor anoxia and increased upper-ocean oxygenation. Our simulated mechanism -a reorganization of global ocean circulation, with reduced importance of northern-sourced waters and a poorer ventilated and deoxygenated deep ocean -has parallels with Pleistocene state transitions in Atlantic meridional overturning (despite a very different continental configuration), and suggests that no simple and predictable relationship between past climate state and oxygenation may exist.

The Late Ordovician mass extinction (LOME) is the first of the 'Big Five' mass extinctions [START_REF] Sepkoski | Phanerozoic marine diversity and the fossil record[END_REF] and was characterized by the disappearance of ~85% of marine species 2 . Extinction clustered in two apparent pulses [START_REF] Harper | End Ordovician extinctions: A coincidence of causes[END_REF] (Fig. 1) -at the beginning of the Hirnantian and near the end of the Hirnantian ('LOME 1' and 'LOME 2', respectively). LOME 1 may be related to habitat loss associated with glacioeustatically-driven sea-level changes or potentially volcanism [START_REF] Harper | End Ordovician extinctions: A coincidence of causes[END_REF][START_REF] Saupe | Extinction intensity during Ordovician and Cenozoic glaciations explained by cooling and palaeogeography[END_REF][START_REF] Hu | Large mass-independent sulphur isotope anomalies link stratospheric volcanism to the Late Ordovician mass extinction[END_REF] , while the drivers responsible for LOME 2 are more debated. In the traditional time-line interpretation, LOME 2 occurred during deglaciation; extinction was then triggered by the expansion of anoxia onto the shelves during the associated sea-level rise [START_REF] Holmden | Nd isotope records of late Ordovician sea-level change-Implications for glaciation frequency and global stratigraphic correlation[END_REF] . However, studies conducted at higher stratigraphic resolution have recently challenged this model [START_REF] Ghienne | A Cenozoic-style scenario for the end-Ordovician glaciation[END_REF][START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF] and suggest that seafloor anoxia [START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF] was already well established during the latest Hirnantian glacial maximum [START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF][START_REF] Stockey | Persistent global marine euxinia in the early Silurian[END_REF] .

The increased areal extent of seafloor anoxia, indicated by δ 238 U, coincident with glacial conditions rather than global warming, presents several apparent conundrums. Firstly, this is seemingly at odds with the well-established conceptual model for typical Mesozoic Ocean Anoxic Events which links volcanism with greenhouse-driven warming and hence the occurrence of global-scale anoxia [START_REF] Jenkyns | Geochemistry of oceanic anoxic events[END_REF] . Indeed, previous ocean biogeochemical modeling studies of Hirnantian climate have reported increased oxygenation during glaciation [START_REF] Pohl | The climatic significance of Late Ordovician-early Silurian black shales[END_REF] -the expected corollary of decreased oxygenation during Mesozoic warming episodes [START_REF] Jenkyns | Geochemistry of oceanic anoxic events[END_REF] . Secondly, despite uranium isotopic evidence for expanded seafloor anoxia globally [START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF] , most preserved stratigraphic sections, deposited at various latitudes in shelf and slope settings, record a shift to more oxic depositional facies during the Hirnantian [START_REF] Melchin | Environmental changes in the Late Ordovician-early Silurian: Review and new insights from black shales and nitrogen isotopes[END_REF] . How can deep ocean deoxygenation develop in a cold climate, and become decoupled from increased shallow marine oxygenation?

We resolve this conundrum through a series of Earth System model experiments that reveal how state transitions in large-scale ocean circulation in response to Ordovician climate change can lead to a dramatic redox restructuring of the marine environment. However, we start by confirming the existence of surface-to-deep redox decoupling during the latest Ordovician and report the results of carbonate I/Ca measured in two Laurentian sections spanning the late Katian to earliest Silurian (Rhuddanian).

I/Ca data and upper-ocean oxygenation during LOME 2

Bulk carbonate I/Ca (Methods) reflects changes in dissolved iodate (IO3 -) concentration, which decreases in the presence of low-oxygen water masses [START_REF] Lu | Iodine to calcium ratios in marine carbonate as a paleo-redox proxy during oceanic anoxic events[END_REF][START_REF] Lu | Oxygen depletion recorded in upper waters of the glacial Southern Ocean[END_REF] . Lower I/Ca values have been documented previously during Mesozoic warming events [START_REF] Lu | Iodine to calcium ratios in marine carbonate as a paleo-redox proxy during oceanic anoxic events[END_REF] , providing evidence for the ocean deoxygenation/warming link. Here, results obtained for two Laurentian sections (Extended Data Fig. 1) show increasing mean I/Ca values from the middle segment of the Hirnantian Isotope Carbon Excursion (mid HICE) to the late HICE and LOME 2 (Fig. 1, Extended Data Fig. 2). The higher I/Ca values recorded in the Laurentian epicontinental seas (water depths < ~200 meters) suggest that latest Ordovician cooling was contemporaneous with increased oxygenation of the upper ocean, which is independently supported by the Anticosti Island organic biomarker analysis of Rohrssen et al. [START_REF] Rohrssen | Lipid biomarkers record fundamental changes in the microbial community structure of tropical seas during the Late Ordovician Hirnantian glaciation[END_REF] and is not easily explained as a diagenetic artifact [START_REF] Hardisty | Perspectives on Proterozoic surface ocean redox from iodine contents in ancient and recent carbonate[END_REF][START_REF] Jones | Sea level, carbonate mineralogy, and early diagenesis controlled δ13C records in Upper Ordovician carbonates[END_REF] (Supplementary Discussion and Extended Data Fig. 3). Furthermore, I/Ca values are higher at Anticosti Island (Canada) [START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF] than at Copenhagen Canyon (USA) [START_REF] Hu | 87Sr/86Sr evidence from the epeiric Martin Ridge Basin for enhanced carbonate weathering during the Hirnantian[END_REF] during LOME 2 (Fig. 1). This is consistent with Copenhagen Canyon being deposited in predominantly poorly oxygenated offshore waters and closer to the large Laurentian oxygen minimum zone (OMZ; Extended Data Fig. 1) core as well as with Fe-speciation data from Vinini Creek (Nevada) recording local anoxic deep waters in an upwelling setting [START_REF] Ahm | Disentangling the record of diagenesis, local redox conditions, and global seawater chemistry during the latest Ordovician glaciation[END_REF] (Fig. 1).

Two models for seafloor anoxia in a cold climate

To reconcile the various observations of increased seafloor anoxia from uranium isotopes [START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF] , but increased upper-ocean oxygenation from I/Ca, we carry out a gridded parameter ensemble of simulations using the cGENIE Earth System model of intermediate complexity [START_REF] Ridgwell | Marine geochemical data assimilation in an efficient Earth system model of global biogeochemical cycling[END_REF] configured with Ordovician continental boundary conditions and reduced (3.65 %) Ordovician solar luminosity value (Methods). We simulate a range of plausible Hirnantian climatic conditions by running the model multiple times using different atmospheric pCO2 values (the first parameter axis in the ensemble), for which we chose a range spanning ×5 through ×24 pre-industrial (280 ppm) concentrations [START_REF] Berner | GEOCARBSULF: A combined model for Phanerozoic atmospheric O2 and CO2[END_REF] -equivalent to 1400 to 6720 ppm.

In the ocean, it is likely that nutrient inventories were lower than modern as a consequence of reduced (compared to modern) Ordovician plant cover [START_REF] Porada | High potential for weathering and climate effects of non-vascular vegetation in the Late Ordovician[END_REF] and high sea level [START_REF] Scotese | PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic (PALEOMAP Project[END_REF] , both tending to suppress terrestrial weathering and cause lower river-nutrient input than today. We therefore also simultaneously considered a range of ocean dissolved phosphate inventories, from 0.2 to 1.0 times the present-day mean oceanic concentration ([PO4]) of 2.159 μmol kg -1 (Fig. 2), which constitutes the 2 nd parameter axis of our ensemble. We hence explore a total of 30 different combinations of pCO2 and ocean phosphate inventory in a regular grid. All simulations were run to steady state (10,000 years) and the ensemble we focus on assumes a reasonable estimate for the Ordovician atmospheric oxygen abundance [START_REF] Lenton | COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time[END_REF][START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF][START_REF] Sperling | Statistical analysis of iron geochemical data suggests limited late Proterozoic oxygenation[END_REF] and hence climatic state (Fig. 2). As expected, the extent of seafloor anoxia (defined as benthic [O2] ≤ 0 μmol kg -1 and the presence of free sulphide (H2S) in the model) shows a generally positive saturating relationship with increasing phosphate inventory -a direct response of the changing bacterial oxygen demand as biological primary productivity and carbon export to the ocean interior (and the strength of the 'biological pump') increase. In contrast, we find an unexpected and highly non-linear relationship between the extent of seafloor anoxia and pCO2, with a prominent maximum in % anoxia reached for intermediate pCO2 values (between ×6 and ×10 CO2) at ×0.4 [PO4].

To identify a plausible scenario through climate-productivity space for the cooling climate leading into the Hirnantian glacial maximum (Fig. 2), we first compare modeled surface ocean temperatures with estimates from carbonate clumped isotopes 26 and d 18 Oapatite 7,27 (Supplementary Results). Hirnantian climate trends are not currently known with sufficient resolution on Anticosti (or elsewhere) to precisely establish the timing of the onset of cooling leading into the late HICE glacial maximum and its relationship to the δ 238 U and I/Ca records. We estimate the magnitude of cooling by comparing paleotemperatures from the Laframboise Member of the Ellis Bay Formation with those from the lower Ellis Bay Formation (because the Laframboise Member is interpreted as representing deglacial transgression rather than the glacial maximum 6 , our estimated temperature changes are best regarded as minimal). Clumped isotopes suggest a decrease from ca. 33 °C to 29 °C, while d 18 Oapatite suggest a decrease from ca. 30 °C to 23 °C.

We use these estimates to posit 2 end-member possibilities for the associated change in atmospheric pCO2 within the model ensemble -a drop from ×24 to ×8.5 CO2 26 (cooling Scenario #1), or a drop from ×7 to ×5 CO2 [START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF][START_REF] Trotter | Did cooling oceans trigger Ordovician biodiversification? Evidence from conodont thermometry[END_REF] (cooling Scenario #2).

Next, and with these 2 plausible cooling scenarios in mind, we seek to constrain pairs of simulations that are characterized by an appropriate expansion of seafloor anoxia in the cold climatic state. Mass balance modeling of the late Hirnantian negative carbonate δ 238 U excursion at Anticosti Island 7 permits us to refine the estimate of anoxia extent by accounting for the marine uranium budget [START_REF] Zhang | Uranium isotopes in marine carbonates as a global ocean paleoredox proxy: A critical review[END_REF] and indicates a shift from near-modern levels of ocean oxygenation (< 1 % seafloor area anoxic) [START_REF] Bartlett | Abrupt global-ocean anoxia during the Late Ordovician-early Silurian detected using uranium isotopes of marine carbonates[END_REF][START_REF] Stockey | Persistent global marine euxinia in the early Silurian[END_REF] , to anywhere from 3 % to 80 % total anoxic seafloor area 8 (Extended Data Fig. 4; Supplementary Results). We can hence further refine our choice of model simulations to pairs in which the extent of seafloor anoxia is initially (in the warm state) close to modern but is substantially expanded in the cold state.

As illustrated in Fig. 2 (and Extended Data Fig. 5), for Scenario #1 (a pCO2 drop from ×24 to ×8.5 CO2), we identify a pair of model experiments consistent with our criteria for anoxia expansion (in this case, from ca. 0 % to 45 % of seafloor area) at a constant nutrient inventory of ×0.4 [PO4]. For #2 (pCO2 decreasing from ×7 to ×5 CO2), an increase in ocean phosphate is required in addition to the cooling in order to trigger deep-ocean deoxygenation. We identify an increase from ×0.2 to ×0.6 [PO4] as consistent with our criteria for seafloor anoxia expansion (Fig. 2, Extended Data Fig. 5). The possibility of an increase in ocean nutrient inventory is supported by models and experiments suggesting that colonization of the continents by (non-) [START_REF] Porada | High potential for weathering and climate effects of non-vascular vegetation in the Late Ordovician[END_REF] vascular 29 land plants could have driven a substantial increase in phosphorus weathering flux during the Hirnantian [START_REF] Porada | High potential for weathering and climate effects of non-vascular vegetation in the Late Ordovician[END_REF][START_REF] Lenton | First plants cooled the Ordovician[END_REF] . Additional simulations run under alternative assumptions of atmospheric pO2 (Extended Data Fig. 6) reveal that these nutrient scenarios are not unique, and that different combinations of pO2 and ocean phosphate inventory may produce similar modelled values of seafloor anoxia. However, the major distinction between the two scenarios remains the same: cooling alone, triggers the expansion of seafloor anoxia in Scenario #1, while cooling in Scenario #2 requires a simultaneous increase in ocean nutrient inventory.

Reconciling shallow-and deep-marine oxygenation trends

To distinguish between the 2 model scenarios, we now turn to evidence for upper-ocean oxygenation changes. A key difference between Scenarios #1 and #2 is the change in the oxygen concentrations simulated in the upper ocean (Fig. 3). In Scenario #1, oxygen concentrations in the upper ocean increase in response to climate cooling because oxygen solubility increases in colder water [START_REF] Pohl | The climatic significance of Late Ordovician-early Silurian black shales[END_REF][START_REF] Garcia | Oxygen solubility in seawater: Better fitting equations[END_REF] while the metabolic rates of marine organisms, and thus primary productivity, decrease [START_REF] Crichton | Calibration of key temperature-dependent ocean microbial processes in the cGENIE. muffin Earth system model[END_REF] . The rate at which particulate organic matter is degraded in the ocean interior also decreases in response to lower temperatures [START_REF] Crichton | Calibration of key temperature-dependent ocean microbial processes in the cGENIE. muffin Earth system model[END_REF][START_REF] John | Temperature-dependent remineralization and carbon cycling in the warm Eocene oceans[END_REF] , with less remineralization (and hence demand in oxygen) occurring in the upper water column, leaving a relatively larger flux to be consumed at greater depths. A more efficient (deeper penetrating) biological pump also reduces the return nutrient flux to the ocean surface, further lowering primary productivity (see Methods and Supplementary Results). In contrast, for Scenario #2, upper-ocean (subsurface) oxygen concentrations decrease primarily as a consequence of a dominant nutrient inventory-driven enhancement of the biological pump which overwhelms the effects of cooling (Supplementary Results). Our I/Ca measurements then distinguish between the two scenarios and identify Scenario #1 as best capturing the overall Hirnantian ocean redox changes, with the spread of seafloor anoxia 7 paired with increased oxygenation of the upper ocean (Fig. 1).

By applying multiple data constraints to an ensemble of different realizations of Ordovician climate and marine biogeochemical cycle states, we can reconcile the divergent trends in seafloor and upper-ocean oxygenation in response to ongoing global cooling. In the model, we account for the observations as a result of a fundamental reorganization occurring in large-scale ocean circulation patterns. In the warmest ×24 CO2 state, deep-water formation in the model predominantly occurs at the highest ocean latitudes (which in the Late Ordovician were located in the Northern Hemisphere) and ventilates much of the Ordovician seafloor with highly oxygenated waters (Fig. 4i-l). As pCO2 is lowered through to ×8.5, seasonal sea ice appears over the North Pole (Fig. 4e-h). Seasonal freshwater fluxes associated with retreating sea-ice cover locally stratify and thus stabilize the water column (Supplementary Results). This mechanism, possibly strengthened by a positive feedback between reduced northwards salinity transport and deep-water formation and hence overturning, leads to a shutdown of deep-water formation over the North Pole. Although Southern Hemisphere overturning strengthens, this is restricted to the West of the low-latitude Laurentia and Baltica continents and leaves much of the ocean floor (particularly to the East of Laurentia and Baltica) poorly ventilated (Fig. 4f-g). We find support for such a shift in large-scale ocean circulation occurring in response to cooling, in the higher resolution MITgcm study of Pohl et al. [START_REF] Pohl | The climatic significance of Late Ordovician-early Silurian black shales[END_REF] . Despite using a different continental reconstruction of the Late Ordovician 34 , they found a shut-down of northern-sourced, deep water in response to a reduction in the solar constant and cooling, with a similar change in the pattern of global overturning (Supplementary Results) and a restriction of southern-sourced oxygenation to the West of Laurentia and Baltica. More idealized basin geometry (and coincidently also using MITgcm) modelling also demonstrates the potential for abrupt transitions in meridional overturning circulation and Rose et al. [START_REF] Rose | The role of oceans and sea ice in abrupt transitions between multiple climate states[END_REF] find this specifically in response to the incipient growth of a polar sea-ice cover. Regardless, further work using coupled climate models that can account for the potential role of inter-annual variability and ocean-atmosphere modes in the statedependence of large-scale circulation patterns in the Ordovician (that we are unable to address with the cGENIE model) is important and would help identify which modes of Late Ordovician ocean circulation are robust (and which might be model-dependent) and how they respond to climate change.

We suggest that a helpful analogy for better visualizing the Ordovician circulation and oxygenation changes that occur in our model between ×24 and ×8.5 pCO2 (Scenario #1), can be found during the glacial-interglacial cycles of the late Pleistocene. Specifically, we note how northern-sourced, deep-water production in the Atlantic decreases as climate cools (and is replaced at depth by deep water sourced from the South) [START_REF] Adkins | The role of deep ocean circulation in setting glacial climates[END_REF] , with the potential for multiple circulation states with abrupt transitions between them [START_REF] Rose | The role of oceans and sea ice in abrupt transitions between multiple climate states[END_REF][START_REF] Ferreira | Linking Glacial-Interglacial States to Multiple Equilibria of Climate[END_REF] . Furthermore, shallowing of the Atlantic meridional overturning circulation coincides with proxy evidence for reduced oxygenation at depth, despite glacial deep-ocean temperatures being colder than during interglacials [START_REF] Jaccard | Large climate-driven changes of oceanic oxygen concentrations during the last deglaciation[END_REF][START_REF] Lu | I/Ca in epifaunal benthic foraminifera: A semi-quantitative proxy for bottom water oxygen in a multi-proxy compilation for glacial ocean deoxygenation[END_REF] .

Finally, an important conclusion we draw from our I/Ca data and simulations is that increased seafloor anoxia is, by itself, an insufficient explanation for the LOME 2. There is no evidence of reduced local upper-ocean oxygenation associated with LOME 2 at Anticosti; if anything there is evidence of increased oxygenation based on I/Ca at this time (Fig. 1). The majority of Anticosti taxonomic last occurrences also occur above the onset of the negative uranium isotope excursion [START_REF] Delabroye | The end-Ordovician glaciation and the Hirnantian Stage: A global review and questions about Late Ordovician event stratigraphy[END_REF][START_REF] Mccracken | Conodont biostratigraphy and paleoecology of the Ellis Bay Formation, Anticosti Island, Quebec, with special reference to Late Ordovician-Early Silurian chronostratigraphy and the systemic boundary[END_REF] , and globally many of the taxa that disappear during LOME 2 occurred primarily in shallow-water habitats that would have remained well oxygenated under almost any reasonable global change scenario. Thus, we suggest that cooling must have interacted with other factors such as changes in nutrient cycling and primary producer communities [START_REF] Melchin | Environmental changes in the Late Ordovician-early Silurian: Review and new insights from black shales and nitrogen isotopes[END_REF][START_REF] Sheets | Graptolite community responses to global climate change and the late ordovician mass extinction[END_REF] , and potentially heavy metal toxicity [START_REF] Vandenbroucke | Metal-induced malformations in early Palaeozoic plankton REFERENCES CITED IN METHODS[END_REF] , to drive the LOME 2 extinction.

CORRESPONDING AUTHORS

Correspondence and requests for materials should be addressed to A. Pohl (alexandre.pohl@ubourgogne.fr) and Z. Lu (zunlilu@syr.edu).

ACKNOWLEDGMENTS

Description of the numerical experiments

We adopted the Hirnantian (445 Ma) continental reconstruction of Scotese and Wright 22

(Extended Data Fig. 1). In the absence of better constraints, a flat bottom ocean floor was used, which is not expected to significantly impact simulated global benthic ocean oxygenation patterns [START_REF] Osen | Sensitivity of Late Permian climate to bathymetric features and implications for the mass extinction[END_REF] (Supplementary Methods). We calculated solar luminosity at 445 Ma (1318.1 W m -2 ) after Gough [START_REF] Gough | Solar interior structure and luminosity variations[END_REF] , and employed a null eccentricity-minimum obliquity orbital configuration, which provides an equal mean annual insolation to both hemispheres with minimum seasonal contrasts.

For each different atmospheric pCO2 assumption, we first ran the FOAM 50 general circulation model experiment for 2 kyrs using these boundary conditions (Supplementary Methods). We then derived the 2D wind speed and wind stress, and 1D zonally-averaged albedo forcing fields required by the cGENIE model, using the 'muffingen' open-source software version v0.9.20, which is assigned a DOI: 10.5281/zenodo.4615664. cGENIE model simulations were initialized with a sea-ice free ocean and homogeneous temperature and salinity in the ocean (5 °C and 34.9 ‰, respectively). The model was then integrated for 10,000 years until equilibrium of the deep-ocean oxygen concentration. Results of the mean of the last simulated year were used for this analysis.

Materials and I/Ca methods

Samples from two Upper Ordovician sections were analyzed for I/Ca: Anticosti Island, Quebec, Canada and Copenhagen Canyon, Nevada, USA (Supplementary Methods). Rock samples were crushed and homogenized to fine powder. Around 4 mg of powdered samples were weighed and

  of 8.38 % -×0.4 the modern value of 20.95 %. Changes in the extent of seafloor anoxia simulated in response to changes in both atmospheric pCO2 and ocean [PO4] are shown in Fig. 2. At ×0.2 [PO4], the model Ordovician seafloor is fully oxic, regardless of the value of pCO2

Figure 1 .

 1 Figure 1. Stratigraphic, glacio-eustatic and geochemical records from (a) western Anticosti

Figure 2 .

 2 Figure 2. Simulated sensitivity of ocean anoxia to climate and marine productivity. Extent of

Figure 3 .

 3 Figure 3. Upper-ocean oxygen concentrations for Scenarios #1 and #2. Upper row (a, c, e)

Figure 4 .

 4 Figure 4. Ocean circulation response to cooling. Meridional overturning circulation (1 st column,

  6, 1-7 (2015). 44. Desrochers, A. et al. A far-field record of the end Ordovician glaciation: The Ellis Bay Formation, Anticosti Island, Eastern Canada. Palaeogeogr. Palaeoclimatol. Palaeoecol. The cycling of carbon and associated tracers in the ocean is based on a single (phosphate) nutrient limitation of biological productivity following Reinhard et al. 47 but adopts the Arrhenius-type temperature-dependent scheme for the remineralization of organic matter exported to the ocean interior of John et al. 33 . Compared to recent Earth System models of the Coupled Model Intercomparison Project phase 6 (CMIP6), the simplified climatic component of cGENIE offers a rapid model integration time, facilitating the numerous sensitivity tests required to quantify the impact of the uncertainty in the model boundary conditions pervading deep-time studies, as well as the 10,000-year-long model simulations necessary to ensure deep-ocean equilibrium (Supplementary Results).

	296, 248-263 (2010).
	45. Finney, S. C., Berry, W. B. N., Cooper, J. D. & Ripperdan, R. L. Late Ordovician mass
	extinction: A new perspective from stratigraphic sections in central Nevada. Geology 27,
	215-218 (1999).
	46. Goldman, D. et al. The Ordovician Period. in Geologic Time Scale 2020 (eds. Gradstein,
	F. M., Ogg, J. G., Schmitz, M. D. & Ogg, G. M.) 631-694 (Elsevier, 2020).
	doi:10.1016/b978-0-12-824360-2.00020-6.
	METHODS
	Description of the model
	Simulations were conducted using the Earth System model of intermediate complexity
	cGENIE 19 (Supplementary Methods). cGENIE is rooted in a 3D ocean circulation model, but for
	speed, includes only a 2D energy-moisture-balance atmospheric component plus sea-ice module.
	The model was configured on a 36×36 equal area grid, with 16 unevenly-spaced vertical levels in
	the ocean.

The authors thank M.J. Melchin, D. Horton and an anonymous reviewer for comments that significantly strengthened the manuscript. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 838373 (to A.P.), from NSF EAR-1349252, OCE-1232620, and OCE-1736542 (to Z.L.), from the Natural Sciences and Engineering Research Council of Canada (NSERC Discovery Grant) and NSF 17-536 (to A.D.), from the David and Lucile Packard Foundation (to S.F.), from NSF grant 1736771 and the Heising-Simons Foundation (to A.R.) and from the National Natural Science Foundation of China (41520104007, 41721002) (to Y.S. and M.L). Calculations were partly performed using HPC resources from DNUM CCUB (Centre de Calcul de l'Université de Bourgogne).

DATA AVAILABILITY

The I/Ca data can be downloaded from Zenodo (https://zenodo.org/record/5136966#.YP5vClMzbu6).

CODE AVAILABILITY

The code for the version of the 'muffin' release of the cGENIE Earth System Model used in this paper, is tagged as v0.9.20, and is assigned a DOI: 10.5281/zenodo.4618023.

AUTHOR CONTRIBUTIONS

A.P., Z.L. S.F. and A.R. designed the study and wrote the manuscript with input from all coauthors. AP. conducted the FOAM and cGENIE experiments and led the analysis of the model results. Z.L., W.L., and R.H. carried out the I/Ca analyses. Z.L. led the analysis of the I/Ca results. R.G.S. conducted the mass balance modeling of the uranium isotope cycle.

COMPETING INTERESTS

The authors declare no competing interests. rinsed with de-ionized water to remove soluble iodine. Nitric acid (3%) was added to dissolve carbonate. The samples were then diluted to a consistent matrix containing ~50 ppm Ca, 0.5% tertiary amine (to stabilize iodine), 5 ppb In and Cs (internal standards). The measurements were performed immediately on a quadrupole ICP-MS (Bruker M90) at Syracuse University. Iodine calibration standards were prepared daily from KIO3 powder. The sensitivity of iodine was tuned to 80-100 kcps for a 1 ppb standard. The reference standard JCp-1 was analyzed repeatedly to monitor long-term accuracy [START_REF] Lu | Iodine to calcium ratios in marine carbonate as a paleo-redox proxy during oceanic anoxic events[END_REF][START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF] . The detection limit of I/Ca is on the order of 0.1 µmol mol -1 .

Configuration files for the specific experiments presented in the paper can be found in the directory: genie-userconfigs/MS/ pohletal.NatGeo.2020. Details on the experiments, plus the command line needed to run each one, are given in the readme.txt file in that directory. All other configuration files and boundary conditions are provided as part of the code release.

A manual detailing code installation, basic model configuration, tutorials covering various aspects of model configuration and experimental design, plus results output and processing, is assigned a DOI: 10.5281/zenodo.4615662.