Carbon, Nitrogen, and Sulfur Elemental Fluxes in the Soil and Exchanges with the Atmosphere in Australian Tropical, Temperate, and Arid Wetlands

: Australian ecosystems, particularly wetlands, are facing new and extreme threats due to climate change, land use, and other human interventions. However, more fundamental knowledge is required to understand how nutrient turnover in wetlands is affected. In this study, we deployed a mechanistic biogeochemical model of carbon (C), nitrogen (N), and sulfur (S) cycles at 0.25 ◦ × 0.25 ◦ spatial resolution across wetlands in Australia. Our modeling was used to assess nutrient inputs to soil, elemental nutrient ﬂuxes across the soil organic and mineral pools, and greenhouse gas (GHG) emissions in different climatic areas. In the decade 2008–2017, we estimated an average annual emission of 5.12 Tg-CH 4 , 90.89 Tg-CO 2 , and 2.34 × 10 − 2 Tg-N 2 O. Temperate wetlands in Australia have three times more N 2 O emissions than tropical wetlands as a result of fertilization, despite similar total area extension. Tasmania wetlands have the highest areal GHG emission rates. C ﬂuxes in soil depend strongly on hydroclimatic factors; they are mainly controlled by anaerobic respiration in temperate and tropical regions and by aerobic respiration in arid regions. In contrast, N and S ﬂuxes are mostly governed by plant uptake regardless of the region and season. The new knowledge from this study may help design conservation and adaptation plans to climate change and better protect the Australian wetland ecosystem.


Introduction
Anthropogenic development, climate change, and land uses such as mining, intense agriculture, cattle farming, gas extraction, and urbanization [1,2] threaten the resilience of wetlands and can increase their rates of degradation.Understanding the impacts of these external forces on greenhouse gas (GHG) emissions and nutrients (carbon C, nitrogen N, and sulfur S) turnover in wetlands is pivotal to develop more effective restoration strategies and future adaptation plans, in particular in Australia, which is considered the second most arid continent [3] and one of the most impacted by climate change [3].For example, extreme climatic events such as cyclones, droughts, intense heat waves, and bush fires followed by floods increase the risk of local vegetation loss.Rising sea level may shift native vegetation towards salt-tolerant mangroves, resulting in a loss of biodiversity and a change in soil nutrient input quality and quantity [4][5][6].Increasing CO 2 concentration in the atmosphere may change wetland productivity, affect biogeochemical equilibria, and alter nutrient storage capacity of wetlands [7].Similarly, changes in wetland hydrology may lead to abrupt releases in C, N, and S in the form of CO 2 , N 2 O, and other gases into the atmosphere, making wetlands a large source of gas emissions rather than a sink of nutrients [8].Agricultural production can also pose nutrient management pressures on Australian wetlands [9].Most agricultural activity is concentrated in the wheatbelt area in the southern and south-western Australia, where seasonal rainfalls alleviate water supply constraints to permit cropping and grazing.From the early 1990s, a transition in land management from crop rotation to fertilization has led to an increase in nutrient load per unit area [10,11], resulting in an increase in N 2 O emissions, N leaching to the groundwater, and production of high-nutrient runoff [12,13].These factors can cause deterioration of the ecosystem services provided by wetlands [14][15][16].
Over the long term, excess nutrient can affect the rate of soil organic matter (SOM) decomposition, while changes in hydroclimatic conditions can amplify CH 4 and CO 2 emissions [17], making wetlands both the trigger and victim in this scenario.Field measurements may help quantify processes and mechanisms affecting a wetland; however, it is unrealistic to implement point-scale distributed measurements over the entire Australian continent.Hence, mechanistic models that integrate the feedback between external forces and biologically-mediated processes have become important tools for understanding nutrient dynamics in soil and improving the quantification of nutrient budgets.Current mechanistic models [18][19][20][21][22][23] integrate hydroclimatic variables, soil properties, land use management, vegetation dynamics, and soil biology to quantify GHG emissions and C sequestration in soil.Despite this, there is a lack of representation of nutrient feedbacks and runoff inputs [24].Specifically, the elemental fluxes of C, N, and S in soil, their bioavailability changes during a year, and how these affect the internal fluxes in terms of emissions, leaching, plant uptake, and stock variability are underrepresented.Hence, a more insightful assessment of the nutrient cycle in Australian wetlands is required to mitigate the effects of human activity and climate change.
The aim of this study is to understand the different allocation of nutrients in each biologically-mediated reaction and shed light on the complex feedback between nutrients and the overall bio-production of GHG.We contend that modelling assessments are the only means currently available to fill in knowledge gaps over geographic scales.Biogeochemical processes are commonly studied through in-situ and laboratory-based experiments, which provide an in-depth understanding of the governing mechanisms.However, they fail to elucidate how the complex feedback between different processes affect long-term soil carbon, nitrogen, and sulfur stocks and emissions, nor can they reveal this critical information at a large spatial scale.To this end, we apply a complex biogeochemical reaction network, BAMS4 (Biotic and Abiotic Model for SOM-version 4) in the BRTSim computational solver to provide a regional scale assessment of C, N, and S cycles in Australian wetlands over a 0.25 • × 0.25 • resolution grid in the year from 2008 to 2017.Agricultural runoff, wet and dry deposition, and nitrogen fixation were coupled to BAMS4 to include the most important nutrient inputs to wetlands.With the BRTSim-BAMS4 framework, we: (1) estimate the spatial distribution of CH 4 , CO 2 , and N 2 O gas emissions from wetlands over Australia; (2) study the nutrient contribution of inputs and their seasonal variability, (3) analyze the elemental partitioning throughout the biogeochemical reactions in soil, and (4) identify the most important processes for C, N, and S cycles.

BAMS4 Biogeochemical Network
The Biotic and Abiotic Model for SOM-version 4 (BAMS4, Figure 1) is a simplification of the carbon-nitrogen-sulfur cycle model (BAMS3) designed in [25], and based on [26][27][28].The C cycle includes the depolymerization of two organic polymer pools (PolyC, PolyCN), aerobic respiration of three organic monomers (MonoC, MonoCN, and MonoCS, representing the organic C, N, and S), CH 4 genesis, and aerobic oxidation.The C cycle was linked to the inorganic N, S, and P cycles to account for nitrification (NH , and C, N, and S protection on the soil matrix (e.g., mineral surface binding).Twelve microbial functional groups were included in the reaction network, including depolymerizing fungi F DEP , heterotrophic prokaryotes B AER , methanogenic prokaryotes B MGB , methanotrophs B MOB , ammonia-oxidizing prokaryotes B AOB , nitrite-oxidizing prokaryotes B NOB , denitrifying prokaryotes B DEN , sulfur-reducing prokaryotes prokaryotes B SrRB , thiosulfate-and sulfide-reducing prokaryotes B ThSRB , sulfate-reducing prokaryotes B SRB , thiosulfate-and sulfide-disproportioning prokaryotes B SDB , and photolithoautotrophoxidizing prokaryotes B SOB .The microbial dynamics, described using Michaelis-Menten-Monod (MMM) kinetics, included C and N immobilization but excluded S and P immobilization because they represent less than 3% of the microbial dry mass.The microbial response is controlled by temperature, pH, water availability, and inhibitory redox conditions such as the presence of O 2 and others oxidants (refer to Tables S1-S6 for relevant inhibition factors and Figure S1).Each pool can be expressed in its aqueous, protected, and gaseous forms.Kinetic reactions were used to describe the SOM protection [26], while equilibrium reactions were used to describe the protection of inorganic species.The reaction network also includes plant uptake of NH + 4 , NO − 3 , SO 2− 4 , and PO 3− 4 and CH 4 aerenchyma transport from the root zone to the atmosphere.
The reaction parameters and corresponding references are summarized in Tables S1-S6 in Supplementary Information.
The MODIS-IGBP land cover product (MODIS, [45]) was used to identify wetlands affected by agricultural runoff (Figure S2a), then the nutrient-rich agricultural runoff to wetlands was calculated using the Runoff Curve Number model (NRCS-CN, [46]).A digital elevation map was used to determine the flow direction (Data Announcement 88-MGG-02, Digital relief of the Surface of the Earth NOAA, National Geophysical Data Center, Boulder, Colorado, 1988).
All data collected for the purpose of modeling are heterogeneous in resolution but were harmonized to the same resolution and bounding box.We implemented two types of resolution harmonization in the MATLAB environment.The first is a linear interpolation for the continuous variables.This was implemented with a conservative procedure when applied to mass or energy flows (e.g., rainfall, solar radiation, NPP, etc.) that consisted of calculating the area integral at the original resolution, and renormalize the maps at the desired resolution by the areal integral so that the total regional flux is conserved.The nearest-neighbor interpolation was used for categorical variables (e.g., Köppen-Geiger climatic regions, etc.), instead, to preserve sharp details.
Table S2 summarizes all the details and references for each database.

Computational Domain and the BRTSim Solver
We selected all the grid cells in the SWAMPS database [30] where the maximum average wetland extension was greater than 5% in the period 2008-2017; hence, the computational domain covers about 8500 grid cells.Wetland area fraction was estimated by coupling remote sensing [29] and inventory dataset from the Global Lakes and Wetland Dataset (GLWD) [47].The wetlands represented by the green grid cells in Figure S2a are the yearly average (from a monthly dataset), and they occupy only a small fraction of the grid cell.The extent of wetland area across the season is well illustrated by Lake Eyre in central Australia, which can vary from almost 0 (dry season) to nearly 10,000 km 2 (wet season).BAMS4 was solved along a 2 m vertical column, which included 3 layers of the vadose zone (30, 30, and 40 cm thickness), a bottom layer of 100 cm, and 4 atmospheric layers above the soil.The general-purpose multiphase and multi-species biorective transport simulator (BRTSim-v4.0e) was used to solve BAMS4 dynamically over space and time.BRTSim includes solvers for dynamic water and heat flow, as well as aqueous and gaseous diffusion-advection, adsorption and desorption in the soil matrix, gas dissolution, and biological dynamics.
The equilibrium soil organic C [33] was used to initialize the SOM soil profile.The organic C was partitioned into each specific pool of SOM following the predicted steady state profile in [26].Each grid cell was run for 250 years re-looping the hydroclimatic variables of the period 2008-2017 to ensure that all the biogeochemical reactions reached a statistical steady state.The final analyses included only the last 10 years of the simulation.

Methods of Analysis
Wetlands were grouped by climatic class (tropical Tr, arid Ar, and temperate Te) using the first level of the Köppen-Geiger (KG) classification (Figure S2b, [48]).We calculated the monthly long-term average in each grid cell and each climatic class.Lastly, we aggregated the months into four seasons, i.e., summer (December, January, February), fall (March, April, May), winter (June, July, August), and spring (September, October, November).
We analyzed the input fractional contribution, the internal fluxes expressed as emissions, leaching (nutrients leaving the root zone), plant uptake, and changes in the size of each SOM pool in each climatic class and season.This analysis was carried out for each of the C, N, and S cycles.
The pools present in the root zone were then analyzed to understand the most abundant among aqueous, protected, solid, and assimilated into the biomass.The gas emissions were analyzed to calculate the annual average GHG emissions from all Australian wetlands for the period 2008 and 2017.
Finally, we analyzed the elemental flux partitioning of C, N, and S throughout the reaction network by tracking the elements in each reaction of BAMS4.

Inputs to Soil
C, N, and S enter the soil via the net primary productivity (NPP), atmospheric deposition, runoff, and N 2 fixation.Here, we present their fractional contribution spatially and over the seasons.
C entering the soil comes from aboveground NPP AG (e.g., lignin and hemicellulose) and belowground NPP BG (e.g., root exudates, mainly simple sugars).The NPP AG is predominant all year and among all the climatic classes, reaching a maximum of 60% in the tropics in fall (Figure 2a).On average, the ratio between NPP AG and NPP BG changes by approximately ±10% across seasons and climatic classes.The N inputs to soil come from NPP AG , NPP BG , atmospheric deposition (Figure S3a), N 2 fixation, and agricultural runoff.The latter contributes only in the arid and temperate regions between 2 and 7% of the total N input.Deposition is the major contributor (almost 70%) of the S input in the tropical and arid regions (Figure S3b), while NPP BG and NPP AG contribute mostly in the temperate region.Agricultural runoff contributes only 2 and 5% of the total S input in arid and temperate regions, respectively.Long term averages for each season and each KG class are reported in Table S3.
The C:N and C:S ratios of the inputs are significantly different across each climatic zone (Table 1).The greatest C:N ratios are found in the tropics associated with tropical forest (on average 15.35), followed by temperate (13.61) and arid (9.51) regions.Even though the average S input in the temperate region is the highest (Figure 3c), the C:S ratio in this area is also the greatest (72.57), followed by tropical (62.09) and arid (37.31) regions.The C:N and C:S ratios vary significantly between the seasons, particularly in the tropical and temperate regions, where they can change by ±25% .

C, N, and S Fluxes in Soil
The SOM undergoes biological respiration producing C gas emissions (CH 4 and CO 2 ), leaching (i.e., nutrients leaving the root zone), or remaining in the root zone and changing the size of the SOM pools, here after called ∆ stock (i.e., aqueous and protected species).The tropical region has the highest C input between fall and winter (dry season), showing the highest change in C stock (light blue bar in Figure 3a).In contrast, the C input decreases in spring and summer, but the emission rate increases, resulting in C stock consumption (i.e., negative ∆ stock).The arid region shows a very low C input (maximum 0.3 g-C m −2 d −1 ).The emission rate remains almost constant across the season, while ∆ stock is negative during summer and fall when the temperature is favorable for biological activity.C input in the temperate region is the highest among the three climatic classes, and no negative stock variation rates occurred.The highest emission rates occur in summer and fall, while they slow down in winter and spring, highlighting the well known hysteresis that links temperature and C emissions [49].Despite a substantial change in gas emission rate and stock variability over the year, the leaching rate of C shows only minor changes in all climatic classes.However, most of the C input to soil undergoes biological respiration, both aerobic and anaerobic, largely contributing to global GHG emissions.
The N input to the soil is almost one order of magnitude less than C in each season (Figure 3b), and most of the N in the soil is removed by plant uptake (green bar Figure 3b).The ∆ stock is always positive, resulting in N accumulation, except during winter in the temperate region.The total N input reaches a minimum in this season, on average 0.1 g-N m −2 d −1 , but gas emission, leaching, and plant uptake rates remain almost constant, degrading the soil N gained in the previous seasons.S input to the soil is almost one order of magnitude less than N (Figure 3c).Similar to N, S plant uptake is the most significant flux followed by leaching, while S gas emissions are almost negligible.Negative ∆ stock occurs during summertime in temperate region only.Fluxes are expressed as grams of element per unit surface area of wetland per day.
Even though Western Australia (WA) has the greatest wetland area (Table 2), Queensland has the greatest CH 4 emissions followed by New South Wales (NSW) and the Northern Territory (NT) due to high input of C as NPP and a favorable temperature.A similar pattern was found for N 2 O and CO 2 emissions.WA also shows high N 2 O emissions resulting from the extra input of fertilizers from agricultural runoff in the wheatbelt (southwest area of WA).Spatially, the east coast, Tasmania (TAS), and the southwest coast of WA represent hot-spots of all the three GHG emissions, releasing annually ≥100 g-CH 4 m −2 , ≥1250 g-CO 2 m −2 , and ≥0.3 g-N 2 O m −2 (Figure 4).The climate in these areas is temperate and characterized by significant rainfall compared to the arid climate of the inland, therefore the wetland area remains constant across the months, ensuring anaerobic conditions.In addition, favorable hydroclimatic conditions enhance vegetation resulting in high nutrients input (Figure 3), resulting in high substrate availability for GHG genesis.Moreover, Tasmania represents the Australian state with the greatest GHG emissions per unit of wetland area.
N 2 O emission hot-spots can be found in high fertilization areas such as the wheatbelt (blue and yellow grid cells in Figure S2a), and where high N deposition occurs such as in TAS and NT (Figure S3a).In contrast, no or limited CH 4 emissions occur in the central part of Australia due to very small wetland fractional area and very low NPP in arid conditions, allowing only short grass or bush to grow [50].However, this area has very high CO 2 emissions (Table 2) because the topsoil is mostly unsaturated during the year, thus allowing aerobic degradation at a rate approximately ten times faster than anaerobic degradation.Tropical and temperate regions have less wetland area than the arid regions, but CH 4 emissions are three and six times higher, respectively.Moreover, with approximately the same wetland area, the temperate region has three times greater N 2 O emission than the tropics due to agricultural fertilization.

Phase Fractions in Soil
The SOM can be present as solid (PolyC and PolyCN), aqueous (dissolved form), protected (bound to the soil matrix), and immobilized in the microbial biomass (only C and N) as represented in Figure 5.
The protected pool of C is the major phase in each climatic class, followed by the solid, aqueous, and immobilized phases.The arid region shows a particular pattern characterized by all the C in the protected phase (Figure 5a).The soil C content is also significantly affected by the climatic condition and land cover.The arid region, which is located in the internal part of Australia and characterized by extreme climatic conditions adverse to plant growth (mostly drylands of low bushes), has limited input of new C to the soil, while wet and cool regions commonly have higher C content due to high biomass production.In fact, the temperate region shows the highest C content, 36 g-C kg −1 soil , while tropical and arid regions have 13 and 17 g-C kg −1 soil , respectively.The tropical region has the lowest C content despite high NPP, because of a high respiration rate.Likewise, N and S are mostly protected or in the aqueous phase (Figure 5b,c).The particulate SOM (PolyCN) and the immobilized N have negligible mass fractions.Despite the significant variability in C:N and C:S ratios in the input quality, the ratio between nutrients in the root zone remains almost constant over the year (Table 1).The temperate region shows the highest average C:N and C:S ratios in the root zone (18.7 and 98.35, respectively), followed by tropical (13.07 and 80.6), and arid regions (10.76 and 46.94,Table 1).

Mass Flux Partitioning
In BAMS4, we tracked and normalized to 100 the mass passing through each individual reaction, including plant uptake of CH 4 , N, and S (Figure 1), to quantify the elemental mass flux throughout the reaction network.
95% of the C mass flux passed through seven processes (Figure 6a), namely depolymerization R1 (between 10 and 20%) and R2 (between 3 and 7%), MonoC aerobic respiration R3 (between 10 and 20%), MonoCN aerobic respiration R4 (between 7 and 10%), MonoC anaerobic respiration R6 (between 30 and 70%), CH 4 aerobic oxidation R7 (between 10 and 15%), and CH 4 aerenchyma transport (between 5 and 3%).The remaining C was used as a substrate for N denitrification and S reduction.The C partitioning mass fluxes in the tropical and temperate regions are very similar.The main difference is that more C passes through PolyC in the tropics.In fact, this class includes forested wetlands, which have a high C:N ratio of the NPP AG and high polymer content (Sections 2.2 and 3.1).In contrast, the arid region shows less methanogenesis (R6) and more aerobic respiration of both MonoC (R3) and MonoCN (R4) than other region because most of the time the soil column is not fully saturated, therefore resulting in little CH 4 emission (Table 2).
Figure 6b shows that N mass flux partitioning changes between the climatic regions and between seasons.The majority of the N mass flux occurred as plant uptake R16 (between 58 and 91%), confirming the importance of this process, as also highlighted in Section 3.2.The tropical region shows a high mass flux through NH + 4 oxidation (R8) in summer (almost 8%) as compared to the other seasons.A significant amount of N then passes through NO 2 denitrification (R11), NO denitrification (R12), and N 2 O denitrification (R13).Denitrification acquires more importance over summer and spring rather than fall and winter, due to the favorable temperature conditions supporting this process.92% of N is taken up by plants via R16 and R17 in arid regions; the remaining 8% passes through denitrification.The N partitioning through BAMS4 in temperate regions is similar to tropical regions, although shifted back by one season (e.g., the partitioning in the summer in tropical regions is similar to that in the spring in temperate regions) due to differences in precipitation and temperature regimes.Tropical and temperate regions show a more heterogeneous distribution in the mass fractions across the reactions than the arid region due to a high soluble N fraction available in the soil, which is immediately accessible by soil microbes (Figure 5b).
The S cycle shows a different behavior depending on the climatic classes and seasons (Figure 6c).All S in arid regions goes through plant uptake (R31), with no changes during the year due to very low aqueous concentration and S soil input (Figure 5c and Figure 3c).In the tropical region, 95% of the S pass through plant uptake, followed by SO 2− 4 reduction (R28), S 2 O 2− 3 oxidation (R23), and HS − oxidation (R24).Between 40 and 85% of the S mass flux occurs as plant uptake in the temperate region, followed by SO 2− 4 reduction (R28), HS − oxidation (R24), S 2 O 2− 3 oxidation (R23), HS − oxidation (R18), and S oxidation (R22).The high variability of reactions that occurs in the temperate region is possible because the S input to soil is on average 0.023 g-S m −2 d −1 (Figure 3c), which is almost two and three times greater than in the tropical and arid regions, respectively.

Discussion
The annual CH 4 emissions quantified in this study is in line with [51], a CH 4 global inter-comparison model assessment that estimated between 0 and 13.1 Tg-CH 4 in Oceania (inclusive of Australia, New Zealand, Fiji, Tonga, French Polynesia, New Caledonia, Papua New Guinea, and Samoa).Australian GHG assessment has been overall neglected because of a paucity of field measurement towers [52].Only few studies report an average emission between 3 µg and 44 mg-CH 4 m −2 h −1 [53][54][55][56][57], depending on the type of wetlands analyzed (swamps, epherameal wetlands, peat, mangroves, saline wetland, forested wetlands, and freshwater wetlands), while with BRTSim-BAMS4 we estimated an average emission of 0.58 mg-CH 4 m −2 h −1 .BRTSim-BAMS4 results were also benchmarked against the longterm annual average land surface temperature data in [58] (Figure S5a), pH, and soil organic carbon (SOC) from the SoilGrids [33] (Figure S5b,c).The results are satisfactory, although pH and SOC are slightly overestimated.
The Australian continent is characterized by three climatic areas, with a substantial variety of rainfall and temperature patterns affecting the dynamic of nutrient cycle and therefore GHG emissions.Temperate and tropical regions have high precipitation rates, leading to high biomass, prolonged flooded time, and high CH 4 emissions.In this area, soil temperature also ranges between 20 to 35 • C between spring and autumn, creating favorable conditions for B MGB prokaryotes to grow [59,60].In the arid region, a lack of precipitation and extreme temperature results in low nutrient input to soil and low GHG emissions.In addition, the surrounding land management and atmospheric deposition can impact the GHG genesis.The excess nutrient can enhance plant biomass (not explicitly accounted for in our model) and microbial activity, in particular, B DEN and B SRB , at the expense of methanogens [61].In fact, we found the highest N 2 O emissions in the wheatbelt region.Moreover, SO 2− 4 is an inhibitor of CH 4 genesis.Its average inhibition factor is 5% across the entire Australia.However, some areas in WA can reach up to 75% inhibition (Figure S4, this factor has to be multiplied by the other inhibitor factors -P = K Ii /(K Ii + C i ), where K Ii is the inhibition constant and C i the concentration).Therefore, the quality of the emission from a wetland is a combination of feedback between microbial dynamics, plant/litter biomass, hydroclimatic conditions, and external nutrients input.
Changes in hydroclimatic conditions and long-term application of fertilizers may alter plant composition [62,63], resulting in major effects on the aerenchyma transport of CH 4 emissions, but also on N and S cycles.A change in plant composition may cause a change in the above-and belowground litter quantity and quality [64], potentially increasing the C input, hence increasing CH 4 or CO 2 depending on the water table [65].Possible drainage of some peri-urban wetlands may decrease CH 4 and N 2 O emissions, but the loss in benefits provided by the ecosystem services to the biodiversity will be enormous.Moreover, future projections forecast an increase in temperature by 0.2 to 1.6 • C [66], which may increase the overall GHG emissions.Sea level rise may affect the northern part of Australia due to high salt intrusion, which may affect the vegetation and microbial activity.In contrast, inland wetlands may be affected by prolonged drought, reducing their area by 20 to 40% [66], hence increasing CO 2 emissions.This loss in areas may also cause ecological problems, in particular for birds nesting and breeding [52].
In this study, we present the first GHG emission assessment in Australia.Although BAMS4 includes many biogeochemical processes, other processes that may affect the emissions are not accounted for, such as salinity and its inhibition on microbial respiration, hence GHG emissions.Future sea level rising and therefore increase of salt concentration should be accounted for in the future modelling development to reduce model output uncertainty.
The scientific community is currently developing different studies to offset GHG emissions, exploiting new agricultural management that can reduce the input of fertilizers and increase the nutrients sinks.The main goal of this research was to analyze in depth the nutrient cycles in wetland soil in Australia, understanding which process overcomes the others, accounting not only for hydroclimatic variables but also anthropogenic input from nutrient-riched agricultural runoff and N and S atmospheric deposition.

Conclusions
We designed a mechanistic biogeochemical model that coupled C, N, and S cycles.Because of the flexibility of the BRTSim-BAMS4 model, we were able to analyze each individual process in each cycle and identify the most important among them.We found that different climatic patterns strongly affect preferential pathways in particular in the carbon cycle.Temperate and tropical regions together produce almost 90% of total CH 4 emission (5.12 Tg-CH 4 year −1 ) throughout Australia, although these come from only 10% of the total wetland area.Instead, arid regions release 51% of the total CO 2 (90.89Tg-CO 2 year −1 ) in the atmosphere due to their dry condition.High temperatures and the specific precipitation patterns (i.e., long droughts that enhance fast aerobic degradation rather than methanogenesis) may increase soil C degradation, even though plant biomass increases.N and S cycles are governed mainly by plant uptake, which removes between 75 and 90% of N and between 90 and 100% of S, throughout Australia.The temperate region has the highest and most heterogeneous input of nutrients to soil and therefore the elemental flow through the C, N, and S cycle is relatively more uniform, and shows the highest N 2 O emissions.N and S fertilization may interfere with the N and S biogeochemistry and, in turn, inhibiting CH 4 emissions by up to 75% in WA.Additionally, N fertilization causes threefold increase in N 2 O emissions in the wheatbelt area, which contributes to an overall emission of 2.34 × 10 −2 Tg-N 2 O year −1 , across Australia.Our study brings more clarity and knowledge of the main biogeochemical processes occurring in wetlands, contributing to information that can conserve and protect unique Australian wetlands that provide many services such as flood attenuation, pollutant sinks, wild-life habitat, and plant biodiversity refuges.As Australia has a great variety of wetlands (e.g., 64 of them are Ramsar-registered and over 900 are nationally relevant ecosystems), we suggest that research effort should be focused on assessing nutrient cycles that underpin wetland productivity and providing tools to create a more effective conservation plan to preserve these ecosystems.

Figure 3 .
Figure 3. Inputs partitioned into fluxes of (a) carbon, (b) nitrogen, and (c) sulfur expressed as gas emission, nutrient leaching, plant uptake, and change in soil pool size (aqueous and protected phases, here labeled as ∆ stock).Tr, Ar, and Te represent tropical, arid, and temperate regions, respectively.The bars (total of all fluxes) represents the C, N, and S input.Fluxes are expressed as grams of element per unit surface area of wetland per day.

Figure 4 .
Figure 4. Annual long-term average emissions per unit of average wetland area between 2008 and 2017 of (a) CH 4 , (b) CO 2 , and (c) N 2 O.

Table 1 .
Figure 2. (a) Carbon, (b) nitrogen, and (c) sulfur fractional input contribution in soil in each climatic class and for each season.Tr, Ar, and Te represent tropical, arid, and temperate regions, respectively.NPP BG , NPP AG , Dep., and Ag.Runoff refer to soil organic matter (SOM) input below and above ground, total deposition, and agricultural runoff, respectively.Long-term C:N and C:S ratios in soil inputs and root zone in different climatic classes and seasons.
C:N of Soil Input (g-C g-N −1 ) C:S of Soil Input (g-C g-S −1 ) C:N in the Root Zone (g-C g-N −1 ) C:S in the Root Zone (g-C g-S −1 )

Table 2 .
Annual long-term average CH 4 , CO 2 , and N 2 O emissions and average maximum wetland area reported for each State/Territory, entire Australia, and per climatic class.