Modeling daily soil salinity dynamics in response to agricultural and environmental changes in coastal Bangladesh

Understanding the dynamics of salt movement in the soil is a prerequisite for devising appropriate management strategies for land productivity of coastal regions, especially low‐lying delta regions, which support many millions of farmers around the world. At present, there are no numerical models able to resolve soil salinity at regional scale and at daily time steps. In this research, we develop a novel holistic approach to simulate soil salinization comprising an emulator‐based soil salt and water balance calculated at daily time steps. The method is demonstrated for the agriculture areas of coastal Bangladesh (∼20,000 km2). This shows that we can reproduce the dynamics of soil salinity under multiple land uses, including rice crops, combined shrimp and rice farming, as well as non‐rice crops. The model also reproduced well the observed spatial soil salinity for the year 2009. Using this approach, we have projected the soil salinity for three different climate ensembles, including relative sea‐level rise for the year 2050. Projected soil salinity changes are significantly smaller than other reported projections. The results suggest that inter‐season weather variability is a key driver of salinization of agriculture soils at coastal Bangladesh.


Introduction
The coastal zone comprises only 3% of the earth's surface and accommodate 60% of the world's population, a figure set to increase to 80% by 2050 [Hyun et al., 2009]. Worldwide, about 600 million people currently inhabit low-elevation coastal zones that will be affected by progressive salinization [Dasgupta et al., 2015]. Soil salinity is a major, and the most persistent, threat to irrigated agriculture many deltas, such as in coastal Bangladesh [Clarke et al., 2015]. In saline soils, crop growth is hampered by salt accumulation in the crop root zone. If the upward salt movement caused by evaporation exceeds the downward gravitational movement, salt will accumulate in the root zone. Salt in the soil interferes with the crop growth when its concentration exceeds the tolerance limits of the crop [Rhoades and Merrill, 1975;Ayers and Westcot, 1985]. Most plants suffer salt damage at salinities equivalent to electrical conductivity of the soil saturation extract (ECe) of 4 dS/m or higher. Plant growth is restricted even though enough water may be present in the root zone, because salt reduces the plant's ability to take up water and can also be toxic to the plant. In Bangladesh, >30% of the net cultivable land is in the coastal area. Of the 1.689 Mha of coastal lands, about 1.056 Mha are adversely affected by varying degrees of soil salinity [Soil Resources Development Institute, 2010]. The agricultural sector of Bangladesh constitutes an important component of the national economy accounting for around 21% of the national Gross Domestic Product (GDP). It also employs a significant proportion (25.7%) of the 35 million people who inhabit the coastal zone. For these people, the growth and sustainability of the agricultural sector is of paramount importance to their own prosperity and survival, as well as the national economy. Even though there is ongoing research and development to create salt-and drought-tolerant crop varieties, expected future sea-level rise, climate variability and extremes such as cyclones and storm surges, pose a real threat to livelihoods and food security in coastal Bangladesh [General Economics Division, 2009]. Lazar et al. [2015] show that climatic change alone is unlikely to reduce agricultural productivity in southern Bangladesh, but did not include soil salinization on their calculations.
A recent attempt to model soil salinity for the entire coastal Bangladesh has been presented by Dasgupta et al. [2015], which developed an empirical multi-regression model using monthly averaged observations and projected soil salinity by 2050 in 69 districts of coastal Bangladesh. This suggests that many areas will have significant increases in soil salinity. While the study by Dasgupta et al. [2015] improves our understanding of salinity changes in Bangladesh's agricultural soils, it is a data-driven approach and therefore constraint by the frequency and spatial distribution of the observation points and assumes no changes on the land uses. Models based on one-dimensional salt and water balances that include changes in land uses exits but are limited to farm size simulations or does not include the mix of crops and shrimps and fish aquaculture, or important sources of salinity such as overtopping. For example, Oosterbaan [2002] computes soil salt and water balances in four seasons in a year at farm scale and does not include the salt flux from flooding or mixed agriculture and aquaculture. The catchment water and salt balance model of Bari and Smettem [2006] and the polder model of de Louw et al. [2011] resolves the daily salt dynamic at regional scale (>1,000 km 2 ) and polder scale, respectively, but are limited to the stream salinity (i.e., salinity at the catchment/polder water surface network) instead of soil salinity at root zone.
Our aim is to develop a quantitative framework which is able to simulate daily salinity changes in the soil, taking account all the important sources and sinks of salt across a low-lying delta (Figure 1), at the scale of O(10 4 km 2 ) including changes agriculture/aquaculture practices and crop types. We demonstrate the proposed framework by applying it to a region of coastal Bangladesh O(19,000 km 2 ). The main contribution of this manuscript is on model integration and model structure validation. Establishing the validity of the structure of any system dynamic model is the first logical validation (i.e., before the accuracy testing) [Barlas, 1996]. We present for the first time the use of emulators to resolve the water and salt balance equations at the root zone and at regional scale. Emulators are approximation models of an experiment or a more complex numerical models or simulators [i.e., Tokmakian et al., 2012]. Emulators are primarily used for their speed and flexibility to represent any system. If the approximation of the system is acceptable, and if the inputs are expected to stay within the limits of the training dataset, emulators can be used for prediction. The proposed model outputs include not only daily soil salinity time series but also the water and salt fluxes from all main drivers. The analysis of which water and salt fluxes are the key drivers on the soil salinity in coastal Bangladesh are outside the scope of this work.
This manuscript describes the integrated soil water and salinity model and its components (Figure 1), including validation against existing soil salinity temporal and spatial observations. First, the study region of coastal Bangladesh is described. Second, the method used to calculate the daily soil moisture deficit and soil salt mass balance is presented in detail: (1) the main equation that this model solves is the salt balance at the root zone but, (2) to solve this equation we need to estimate the salt and water flux from ground water due to capillary rise, (3) the inter-dependent daily water balance at soil surface and (4) root zone, and (5) the salinity of the standing water. To solve this set of equations across coastal Bangladesh, we need a good representation of ground water depth and ground water salinity and water and salt fluxes from river and coastal flooding. Process-based models such as Delft3D [WL Delft Hydraulics, 2007], Finite-Volume Community Ocean Model (FVCOM) [Chen et al., 2003], and SEAWAT-MODFLOW [Langevin et al., 2003] provide a good representation at the scale of this study but are too demanding computationally to be run for decadal and longer time scales. Hence, we explore the use of emulators as a way to simulate efficiently the groundwater and estuarine hydrodynamics. The proposed method is then validated against temporal and spatial soil salinity observations. Using the proposed model, soil salinity is projected to 2050 for different climate scenarios for the south-west coastal zone of Bangladesh and compared with soil salinity in year 2009. This area is low-lying: the land elevation ranges from 1 to 3 m above sea level, but most of the study area is poldered (i.e., enclosed by embankments to protect the area from floods and salinity intrusion). The tidal range varies between 0.5 and 4.5 m. The land cover in 2010 was dominated by agriculture (45%), followed by natural vegetation (12%), aquaculture (11%), water (8%), and wetland (8%). This deltaic region provides a range of important ecosystem services, which make it highly suitable for agriculture which provides livelihoods for most of the coastal population. The delta plain of the Ganges-Brahmaputra-Meghna river system supports numerous ecosystem services and livelihoods [Nicholls et al., 2016]. About 85% of the people of the coastal zone are partially or wholly engaged in agriculture. Because of the high population density, over 50% of the households are practically landless having less than 0.2 ha of land. Fishing, crop agriculture, shrimp farming and tourism are the area's main economic activities.
There are three distinct seasons in Bangladesh agriculture: Rabi season (November-March; cool, dry winter), Kharif-1 season (March-June; hot humid summer), and Kharif-2 season (June-November; monsoon). In the 1990s, farms practiced mono cropping (i.e., having only one season crop), but wherever irrigation water is available, crops have been recently started to be cultivated in two and three cycles per year. Multi-cropping has become more common because of a higher awareness of alternative crops due to agricultural extension services and non-governmental organisation (NGO) interventions. Additionally, a higher availability of irrigation water, through the installation of diesel-driven tube-wells funded by NGO loans has enabled the capital investment required for high yielding varieties (HYVs). The traditional crop is rice (varieties are Boro, Aus, Aman, cultivated generally in the Rabi, Kharif-1, and Kharif-2 seasons, respectively), but cash-crop production is increasing and includes crops such as wheat, chilli, potato, mustard, tomato, and grass pea. Furthermore, because of agricultural research and development projects, the traditional local varieties of crops have been almost completely replaced by more resilient hybrid and HYVs.
The highest average maximum temperature in the study area is 33 ∘ C and above during March-May, and the lowest average minimum temperature is about 15 ∘ C in December and January. The south-western region of Bangladesh receives an average rainfall of about 1,730 mm per annum, of which about 78% falls within the 4 months of the monsoon. Monsoon rains are important for both providing soil moisture and irrigation water, and flushing the salinity from the soils. However, soil salinity is more spatially variable due to localized environmental processes and management practices (e.g., irrigation, polderisation, etc.). Precipitation generally increases from west (∼1,700 mm/year) to east (∼2,600 mm/year) in the study area, although there is smaller north (2,300 mm/year) to south (2,600 mm/year) gradients. Figure 3 shows the flow chart of the proposed framework. At the start of the simulation, the user select the inputs; climate ensemble, socio-economic scenario, start, and end dates of the simulation. The climate inputs contains daily information of rainfall and temperature for the entire region and climate-coherent upstream river water flows and sea level time series along the coastline of the study site. The socio-economic inputs contains, among other, information of different cropping pattern types and area extent for different time epochs. The start and end dates of the climate and socio-economic inputs constrains the start and end dates that the model will be able to simulate. Once the user has selected the inputs, these are loaded as boundary conditions for the simulation and used to emulate the delta hydrology. The first emulated variable is the river water elevation at a number of virtual river stations using upstream river discharge and sea level time series as predictors ( Figure 1 shows the location of all 105 virtual stations used). Because the outputs of this soil salinity model are given at union level, a union-specific river elevation is estimated using the virtual river stations; for each union, all river stations within a 50 km distance from union centroid are weight averaged using weights that are inversely proportional to their squared distances from the virtual stations. This weighted union-specific river elevation will be later on used to modulate the surface drainage daily rate and the overtopping flood rate. The upstream river water discharge for the three tributaries and the statistics (mean, maximum, standard deviation) of daily sea level was used as predictors to emulate the river salinity at the same virtual river stations used for river elevation. A union-specific river salinity is calculated using the weighted averaging method described above. The union-specific river salinity is used to calculate the flux of salt due to overtopping and irrigation if river is used as an irrigation source (irrigation source is defined as a cropping pattern input  referred as 1st layer) and salinity at ∼100 m (hereinafter referred as irrigation supply layer) depth at each of the 653 unions are emulated. The predictors used to emulate the depth to groundwater are the pumping rate, sea elevation, and river elevation. Depth to groundwater, together with the calculated water logging and root depth are used to estimate the water and salt flows due to capillary rise. The predictors used to emulate the salinity at different depths were the river elevation at the virtual river stations. Once the main hydrological water and salt fluxes are emulated, together with the loaded climate boundary conditions, the salt and water balance at different reservoirs are calculated. First, the water and salt balance at the soil surface (Equations (3) and (6)), then the water balance at the root zone (Equation (5)) and finally the salt balance at the root zone (Equation (1)). Daily time series of soil salinity at the root zone for each crop, within each cropping pattern are calculated following this flow chart twice; for protected and non-protected area within a union (see Table S4, Supporting Information). The only difference on the drivers of soil salinity between the protected and non-protected area is the way that daily flood rate is estimated. For protected land and non-protected land, flood rate is obtained from emulated inundation depth simulations trained with a hydrological simulation model. In addition to daily soil salinity time series, the model also output the daily water and salt fluxes that drives soil salinity. A union-specific daily soil salinity time series are obtained by area averaging all cropping patterns at each union.

Proposed Modeling Framework and Soil Salinity Governing Equations 2.2.1. Proposed Framework for Soil Salinity Calculations at Regional Scale
In this work, we have used three different climate ensemble named Q0, Q8, and Q16. Daily precipitation and temperature values for 1981-2098 were provided by each ensemble of simulations of the HadRM3P (PRECIS) regional climate model developed by the Met Office Hadley Centre . The model is based on the HadCM3 global climate model and dynamically downscaled to assess regional climate variability. The model runs at a 0.22 ∘ × 0.22 ∘ resolution (∼25 km by 25 km), with 19 vertical levels and 4 soil levels. The ensemble member (labeled Q0) represents an unperturbed version of the model physics. The projections are driven by the Special Report on Emissions Scenario (SRES) A1B scenario, which lies between RCP6.0 and RCP8.5 in terms of atmospheric CO 2 concentrations and global temperature projections. This represents a future for the coastal areas in Bangladesh, which is approximately 4 ∘ C warmer and 9% wetter by the year 2098. This paper considers results from the unperturbed model Q0 and two alternate scenarios to examine the sensitivity of the models developed: one which is warmer and drier (Q8) and another which is warmer and wetter (Q16) over the model domain. To eliminate extreme daily precipitation values in HadRM3P dataset, we represent this spatial distribution with three climatic regions (NW, NE, and SE) by averaging the 25 km × 25 km grid values. NW is Khulna division, NE is the northern part of Barisal division, and SE is the southern part of Barisal division. The upstream water flows for the three main tributaries (Ganges, Brahmaputra, and Meghna) are computed with the Integrated Catchment Model (INCA) model [Whitehead et al., 2015]. The tidal water levels at the sea boundary is provided by the simulation of the Bay of Bengal by Kay et al. [2015] using the Global Coastal Ocean Modelling System (GCOMS). Both, upstream water flows and Bay of Bengal time series, coherent with the Q0, Q8, and Q16 climate ensembles have been used.
Evaporation under standard conditions from the land surface (ET c,t ) is calculated as the product of the Potential evapotranspiration (PE) and the crop coefficient (K c ) , which varies with the environmental conditions, crop type, and development stage. The crop evapotranspiration under standard conditions is the evapotranspiration from disease-free, well-fertilized crops, grown in large fields, under optimum soil water conditions, and achieving full production under the given climatic conditions [Allen et al., 1998]. The crop coefficient (K c ) represent the ratio of (ET c,t ) and PE. For a given crop, (K c ) varies with the stage of development. For bare soils, we assume that K c = 0.8, which is the lowest value recommended by Doorenbos and Pruitt [1977] to be used in tropical regions. The PE was estimated using the same approach followed by Clarke et al. [2015]. Using the FAO Penman Montieth approach with historical average climate data from the CLIMWAT 2 database [Muñoz and Grieser, 2006], daily values of PE were obtained by linear interpolation between monthly values. Inter-annual variability of PE is known to be small and only has a small impact on vegetation development [Fatichi and Ivanov, 2014], so the historical data were assumed to remain constant over a 20-year time period 1981-2000. Estimates of future PE were generated by perturbing the historical climate data for the years 2014, 2030, 2050, and 2080[Clarke et al., 2015. This approach was used in place of calculating daily PE values from HadRM3P-generated climate data because the HadRM3P daily temperature and humidity sequences were generated independently and produced inconsistent sequences of PE. The main impact of climatic change is to raise the mean monthly evapotranspiration rate by approximately 0.15-0.2 mm per day by the 2080s. Annual total PE rises in this period approximately 6% from a baseline value of 1,193-1,263 mm.
We have used the cropping patterns coherent with the socio-economic scenario "Business as Usual" of BAU [Allan and Barbour, 2016]. The BAU scenario is defined as the situation that might occur if existing policies continue and development trajectories proceed along similar lines to the previous 30 years. The model includes a library of 94 different crops (25 species of fish and shrimps has been also included as crop types) organized in five different cropping patterns for each Upazilla (i.e., sub-district). A review of observed cropping patterns (i.e., sequence of crops on agriculture fields) were carried out to identify typical crops, their varieties, and use in coastal Bangladesh. This work utilized the Soil Survey Reports of Bangladesh and the five most frequent cropping patterns were selected for each Upazilla. It is assumed that the cropping patterns are the same in each union within a specific Upazilla. The percent area used for each cropping pattern is as observed. If more than one crop was grown in a cropping pattern in the same season, their percentage area is assumed to be equal. Observed cropping patterns of the early 1990s are used for simulation years before the year 2000, whereas observations of 2009-2010 are used as present day cropping patterns from 2000 to 2014. It is assumed that the cropping patterns are going to be the same as the present day in the future. The properties of the observed crops (rooting depth, crop coefficient, evaporation depletion factor, salinity tolerance, etc.) are partially collated from field observations (Bangladesh Agriculture Research Institute, BARI, Bangladesh Rice Research Institute, BRRI, and Department of Agricultural Extension, DAE, datasets), partially based on a model calibration exercise described in Lazar et al. [2015]. To anticipate future crop varieties, properties of future crops (potential yield and salinity tolerance) were modified based on information published in "Agricultural Technology for Southern Region" report [Bangladesh Agricultural Research Council, 2013]. However, basic crop properties that affect water uptake and other tolerances (e.g., temperature) in the model were not changed, thus these are irrelevant to this paper.
In this study, we used the partial least square (PLS) regression to build a set of linear emulators of river water elevation and salinity at the 105 virtual river locations and ground water depth and salinity at each one of the 653 unions. PLSs regression is a technique that combines the principal component analysis (PCA) with the multiple linear regression [Thompson, 2005]. Its goal is to predict a set of dependent variables from a set of independent variables or predictors. This prediction is achieved by extracting from the predictors a set of orthogonal factors called latent variables, which have the best predictive power. PLS regression is particularly useful when we need to predict a set of dependent variables from a (very) large set of independent variables (i.e., predictors). Since the dimension of the emulators outputs for this study are large, being 105 for the river elevation and salinity emulators and 653 (i.e., one for each union) for ground water depth and salinity and inundation depth emulators, we used the PCA as a pre-processing technique to reduce the number of output dimensions (see Text S1 for a more detail description on emulators).
The water river elevation and lateral water inflow due to river flooding is obtained from an emulator trained with simulated flooding with Delft3D for a set of 14 climate scenarios and sea level scenarios (see Table S2 for details on the simulated scenarios). Delft3D-flow is a hydrodynamic model, which describes temporal variation of flow parameters in three dimensions. The model solves the unsteady mass and momentum conservation equations by using a finite difference method [Hydraulics, 2011] and computes water elevation, water depth, velocity, and discharge at each schematized computational grid point (including land, river, and sea) in the model domain. The domain is bounded in the north by the major rivers of the system, i.e., the Ganges, Brahmaputra, and Meghna. A set of nested regular grids have been used to run Delft3D. The spacing of the grid is on the other side (1 km) in the open ocean and it is reduced to spacing of order (100 m) in the land side. Transition from coarse grid to fine grid is kept smooth to ensure that it does not create any numerical diffusion. For further details about this application of Deltf3D to the study area, see Haque and Rahman [2016]. Main Delft3D outputs are time series, at 10 min time interval, of river elevation at 105 locations and gridded inundation water depth. The Delft3D gridded outputs are aggregated into two different time series for the protected and non-protected land within each union. For each union, the daily inundation water depth (mm) and inundation area (ha) are obtained. The lateral inflow due to flooding in mm/d is obtained as the elevation difference between the day of interest and the previous day. If the elevation difference is positive-inundation depth has increased-the net lateral influx is equal to the elevation difference divided by the time interval. If the elevation difference is negative-inundation depth has decreased-there is no net lateral inflow for this day.
The river salinity is obtained from emulated daily river salinities at the same 105 locations where river elevation is emulated. As for river elevation, a Euclidian distance average value is assigned for each union. The river salinity time series at each one of the 105 virtual stations are calculated using an emulator trained with simulated salinities using the FVCOM [Chen et al., 2003[Chen et al., , 2006. FVCOM is an unstructured-grid, finite-volume, free-surface, 3D primitive equation coastal ocean circulation model. The triangular grid resolution varies from large elements of order (10 km) in the open ocean (maximum 26 km), and small elements of order (100 m) within the narrow river channels (minimum 47 m). It solves momentum, continuity, temperature, salinity, and density equations on a triangular grid using terrain following co-ordinates in the vertical. FVCOM was implemented over the delta extending from beyond the shelf-break in the Bay of Bengal at the south, to the limit of tidal movement in the north. The model is forced with hourly water levels and daily temperature and salinity at the ocean boundary provided by the GCOMS model [Kay et al., 2015]. River discharge is applied as an upstream boundary; a volume of freshwater is introduced to the model at a single grid point at the northern boundary of the model. The daily discharge volume rate comes from INCA [Whitehead et al., 2015]. The tidal and fluvial hydrodynamics were validated against water level observations and tidal constituent analysis [Bricheno et al., 2016]. River salinity was calibrated against observations made by Jahan et al. [2015]. See Table S3 for details on the FVCOM-simulated scenarios. The depth to ground water and the salinity of groundwater are obtained from an emulator trained with MODFLOW-SEAWAT simulations for the period 1981-2098 under three different climate scenarios. MODFLOW-SEAWAT is a software that simulates three-dimensional variable-density groundwater flow coupled with multi-species solute and heat transport [Langevin et al., 2003;Harbaugh, 2005]. The SEAWAT model was set-up for a domain larger than the study region encompassing the greater southwest region in Bangladesh. The model uses groundwater pumping volume, general head, flow and river boundary conditions and includes all the significant rivers in the southwest coastal zone to allow for river-aquifer interaction. The model is conceptualized (and subsequently discretized) as consisting of the soil layer, the first (upper) aquifer (a composite aquifer consisting of complex patterns of very fine to medium sand), an aquitard and the second (deeper) aquifer. The model is discretized into a 2 km × 2 km grid, with 10 vertical layers in the first aquifer and 8 vertical layers in the second aquifer. This was done based on lithological analysis of nearly 2,700 bore logs for the Ecosystem Services for Poverty Alleviation Deltas project [Nicholls et al., 2016]. The main model outputs are the depth to groundwater tables (or piezometric heads) and salinities in the different layers. The SEAWAT gridded outputs are aggregated into an area averaged union time series. The simulated salinity of the first aquifer (less than 10 m depth below water surface) is used to characterize the ground water salinity (C g ).

Main Governing Equations of Soil Salinity
The salt balance at the root zone, which extends from the land surface to the crop-specific maximum root depth, is a statement of the law of conservation of mass. In its simplest non-stationary form, it is: where y rz,t is the salt mass content in (gr) at the root zone at a given time (t) and y rz,t−1 is the salt mass content at the previous time step; IC i is the net salt flux due to infiltration from land surface (I being the surface water infiltration rate (mm/day) and C i being the salt concentration of infiltrated water (gr/l)); GC g is the net salt flux due to capillary rise from ground water (G being the water capillary rise rate [mm/day] and C g being the salt concentration of ground water [gr/l]); IrrC irr is the net salt flux due to irrigation, being Irr is the irrigation rate (mm/day) and C irr is the salt concentration of irrigation water (gr/l); and RC rz is the net salt out flux due to percolation, being R the water percolation rate (mm/day) and C rz is the salt concentration at the root zone (gr/l). To obtain this simple salt balance equation, we make the following assumptions: 1. All salts are highly soluble and do not precipitate; 2. The amount of salts supplied by rainfall is negligible; 3. The amount of salts supplied by fertilizers and exported by crops are negligible.
We can regard the amount of salt at any given time in the root zone (y rz,t ) as being dissolved in the soil water. Because the downward movements of water and salt generally take place at water contents near field capacity, we can logically consider y t to be dissolved in an amount of water W rz,fc , which is the amount of soil water at field capacity expressed in mm. W rz,fc can be determined from: where is the soil effective porosity (typically 0.12 for the silty clay soils of coastal Bangladesh) and D is the depth of the root zone (mm). For the sake of simplicity, the depth of the root zone is defined for each crop PAYO ET AL. SOIL SALINITY IN COASTAL BANGLADESH type as the maximum depth at the end of the development stage. The salt concentration in the root zone is then defined as C rz = y rz,t /W rz,fc and Equation (1) can be now iteratively solved if all the salt and water fluxes ( Figure 1) and the initial soil salt content at the root zone are known. To avoid negative soil salinities when solving Equation (1) the salt removed by percolation at each time step is obtained as the minimum between the amount of salt before percolation (y t−1 + IC i + CRC g ) and the salt carrying capacity (RC r )(RC r ). If the salt carrying capacity due to percolation is smaller than the amount of salt in the soil some will remain and vice versa.
In irrigated areas, during intervals in irrigation or during fallow periods when there is no downward flow or percolation water, water can move upward by capillary forces. The water will be taken up by the roots or evaporate at the surface and salt will accumulate in the root zone or in the soil top layer. The capillary upward flux varies with soil type, depth of water table, and soil water gradient [Ritzema, 1994]. Most top soil in coastal Bangladesh is made of silty clay soils [Clarke et al., 2015]. Because very detailed experiments will be required to determine the groundwater contribution under field conditions, we used the capillary rise estimates for clay soil types suggested by Doorenbos and Pruitt [1977]. No capillary flow occurs if the groundwater table is within the root zone or waterlogging occurs (i.e., standing water depth at soil surface is larger than 0) and the plants cannot form air ducts. The time step of one day implies that some of the calculated flow rates may result in moisture contents exceeding the total water storage capacity at the root zone, limiting the values of capillary rise rate suggested by Doorenbos and Pruitt [1977]. Changes in depth to groundwater due to capillary rise are assumed to be negligible relative to seepage from rivers and groundwater pumping. (1) is the recharge into the root zone, its value is related to the inflow and outflow components of the surface water balance. These components are:

Because the rate of infiltration (I) in Equation
1. Water that reaches the land surface from precipitation; 2. Water that enters the water balance area by lateral surface inflow and leaves it by lateral surface outflow; 3. Water that evaporates from the land surface.
Infiltration in the root zone can therefore be expressed by Equation (3): where I t is the infiltration rate for the time interval Δt (mm/d), it varies between 0 and a union specific maximum infiltration rate, P t is the precipitation rate (mm/d), ET c,t is the evaporation from the land surface under standard conditions (mm/d), F t is the lateral inflow due to flooding (mm/d), DR s,t is the lateral outflow due to surface drainage (mm/d), it varies between 0 and a union and time specific maximum rate, ΔW s,t is the the change in surface water storage (mm). In Equation (3) we have assumed that irrigation only occurs when there is no standing water on the field (i.e. all irrigated water flux percolates to the root zone).
The daily time step implies that the calculated infiltration rates may exceed the soil maximum infiltration rate. If maximum soil infiltration rate at any given time step exceeds the maximum rate, the infiltration rate is assumed equal to the maximum rate and the excess is added to the surface water storage. If there is water standing at the soil surface-water ponding-and the calculated infiltration rate from Equation (3) is less than maximum soil infiltration rate, the infiltration rate is increased until either the maximum is reached, or there is no more standing water at soil surface for the given time interval. The surface run-off rate for each time step is assumed to be equal to the minimum of the maximum surface drainage rate and the excess net inflow of water. If there is still a net water inflow after subtracting the infiltration and evapotranspiration out flux, the excess is assumed to be lost as surface run-off until the maximum surface run-off rate is reached. The standing water level will rise if the net water influx (precipitation, flooding) are larger than the net water out flux (infiltration, run-off, evapotranspiration). Infiltration calculations under pond aquaculture are treated slightly different. For aquaculture, ponding is assumed to be maintained while the aquaculture is active. During the pond active period, infiltration is always maximum, pond salinity is assumed equal to the optimal crop growing salinity and run-off is assumed zero. When aquaculture is finished, run-off is non-zero and calculated as for any other crop.

10.1002/2016EF000530
Most surface drainage systems in coastal Bangladesh are gravity driven and therefore the maximum surface drainage rate depends on the elevation difference between the union ground elevation and the nearby river elevation. The union hypsometric curves (area at a given elevation) are used to assess whether there is a reduction of drainage or not. For each time step, the maximum drainage rate is calculated as: where DR smax,t is the union-specific maximum surface drainage rate (mm/d) at time t and f R is a drainage correction factor [0, 1]. The drainage correction factor (f R ) is assumed to be proportional to the area below the river elevation. If none of the union area is below the nearby river elevation level, the drainage correction factor is 1 and the maximum drainage rate for this time interval is equal to R smax . If all union area is below the nearby river water elevation, the drainage correction factor is 0 and no run-off occurs for this time interval. The river water elevation is defined as the weighted averaged river elevation for each union using the emulated river elevation at 105 locations within the study zone.
If the root zone soil moisture content is above field capacity, the excess water percolates deeper into the subsoil. The percolation rate from the rooted zone can be calculated as: where W rz,t is the soil moisture amount in the root zone at a time step t (mm) and W rz,fc is the soil moisture at root zone at field capacity (mm), and G and I are the capillary rise rate from groundwater (mm/d) and infiltration rate from soil surface (mm/d).
Not all crops are irrigated, but for those irrigated crops, the irrigation needs (mm/day) are calculated as the difference between the readily available water (RAW) and the moisture content on each day. Where the soil is sufficiently wet, the soil supplies water fast enough to meet the atmospheric demand of the crop, and water uptake equals ET c,t . As the soil water content decreases, water becomes more strongly bound to the soil matrix and is more difficult to extract. When the soil water content drops below a threshold value, soil water can no longer be transported quickly enough towards the roots to respond to the transpiration demand and the crop begins to experience stress. The fraction of the total available water (TAW) that a crop can extract from the root zone without suffering water stress is the readily available soil water. TAW is the amount of water that a crop can extract from its root zone, and its magnitude depends on the type of soil (effective porosity) and the rooting depth (varies with crop type and development stage). The depletion fraction is assumed to be only a function of the crop type (i.e., it does not vary with the evaporation power of the atmosphere). At every time step, the soil moisture is firstly calculated assuming no irrigation is needed. For each crop on a given day, if irrigation other than rainfall is needed, and the resulting soil moisture is less than the threshold (we use the RAW threshold), then the amount of water needed to raise the soil moisture level up to the field capacity is calculated. The amount of water needed to raise the soil moisture is the irrigation rate used in Equation (1). The percolation rate is then calculated again by Equation (5) but using the modified soil moisture if irrigation has occurred.
Different irrigation sources have different salinities, which has been found as important factor by Clarke et al. [2015] to assess the soil salinity evolution in coastal Bangladesh. When a crop is defined as irrigated, the irrigation source can be selected among rainfall, shallow ground water, deep ground water, river water or an average of shallow groundwater and river water salinities. If rainfall is selected as irrigation source, the salinity of the irrigation water is assumed negligible in Equation (1) and is not calculated as explained above but is an input. If ground water is selected as the preferred irrigation source, this can be from the first aquifer (∼<10 m depth) or from the second aquifer (∼70-120 m depth) from the river, or from a combination of the river and first aquifer.
The salinity of the standing water (C i ) is defined as y s,t /W s,t the ratio between the salt mass content at soil surface (y s,t ) (y s,t ) and the standing water level (W s,t ). The amount of salt at soil surface at each time step t is calculated as: where RivSal is the salinity (gr/l or ppt) at the union-nearby river stations.

Linear Emulator Accuracy
The emulators' prediction accuracy is evaluated using the root mean-square error (RMSE). To assess the robustness of the emulators, different percentages of available simulation data set are used to train the emulators and the remaining data set is use to validate the prediction accuracy. The simulation points used to train the model are randomly selected and the RMSE is calculated for the data set not used for training. This is repeated 30 times to obtain a mean RMSE for each percentage of data used for training. In this section, as an exemplar, only the accuracy of mean river water elevation is presented. A detail description of the scenarios used to train the emulators and the accuracy of the linear emulators can be found in the Supporting Information (S1 and S2).
The predictors used to emulate the river elevation are the INCA-simulated upstream river discharge rate for each one of the three main rivers, and the maximum, mean, and standard deviation of the daily sea elevation [Whitehead et al., 2015]. Figure 4 illustrates the daily variability of inputs used for the 15 scenarios altogether. The daily statistics (maximum, mean, and standard deviation) of sea elevation were used to capture the hourly variability that is embedded into the Delft3D/FVCOM daily aggregated results. Most of variability, 99.5%, of the river elevation can be explained with the first 30 principal components for all the 105 river locations reducing the emulator dimension from 105 to 30.
Linear emulators can accurately reproduce the simulated mean river elevation ( Figure 5). The mean RMSE when only 1% of the data is used is less than 0.4 m and this error stabilizes around 0.35 m when more data are used for training. A detailed investigation of the emulator behavior reveals that the scenario #14 (i.e., the most extreme temperature and sea level scenario) is the one that consistently produces the largest RMSE. Emulated mean river elevation values are well within a factor 10 uncertainty around the target values and are significantly smaller for larger values of water elevation (>∼2 m) than for smaller river elevation values.

Building Confidence on Soil Salinity Model Temporal Variability
To build confidence on the proposed soil salinity model reproducing the temporal variability of the soil salinity, we have compared the model results against the soil salinity observations reported by [Rahman et al., 2013] under Aman rice and shrimp farming at village level in Satkhira district and the observations of Mondal et al. [2001] at two villages in Khulna district. Both studies reported soil salinity and some components of the drivers of soil salinity included in the simulation model but not all. A comparison between observed temperature and rainfall at Satkhira and Khulna stations against the weather ensembles shows that ensemble Q8 is the closest to the observations, but the daily temperature for the climate ensembles used in soil salinity calculations tends to systematically over-estimate the maximum monthly observed temperatures and under-estimate the lowest temperatures ( Figure S2). Rainfall is well within the observed range with some year to year and inter-climate ensembles variability. Because of these differences between model inputs and observed weather conditions, model confidence is assessed against model behavior relative to the drivers for Q8 ensemble instead that based only on good statistical agreement with observed soil salinities.
The proposed soil salinity model can reproduce the inter-annual variability reported by Rahman et al. [2013], but it is very sensitive to the assumed area coverage dedicated to shrimp and rice. They reported the annual maximum and minimum soil salinity values for period 2010-2012 at Shuktia village level under integrated Aman rice and Shrimp farming-traditionally known as Gher. Gher is a modified rice field with high dikes to keep water inside to cultivate shrimps/prawns. This village is within the Khalishkhali union of 38 km 2 of total area, from which 98% is protected land. The field for cultivating rice in the winter season relies on groundwater irrigation, whereas shrimp are raised in the summer and rainy season using saline water from the nearby Dolua River. Aman rice is also cultivated from July to August. The rice variety is not specified thus we have assumed that a high yield variety of Aman rice is grown from July 1st and last 165 days (including nursery and planting) and Bagda is raised on ponds of 0.55 m depth and with a water salinity of 12 ppt starting in November 15th for a period of 125 days. The calculated soil salinity for Aman rice shows the expected pattern of soil salinity increase during the dry season followed by salt decrease during the monsoon season ( Figure 6). During the rainy season, simulated values are consistently lower than reported annual minimum values.
Step salinity increases during the Bagda raising period are due to pond recharges to maintain the desired pond water depth. Because Rahman et al. [2013] did not reported the area coverage dedicated to shrimp and rice, an assumption must be made to be able to compare the model with the observed soil salinity. Two area averaged soil salinity values are calculated: one value assuming shrimp and rice area are equal (i.e., 50% weight for each crop) and other assuming rice area is slightly larger than shrimp area (70% rice and 30% shrimp). The area-averaged values during the raising of the crops are in good agreement with the maximum and minimum observed annual soil salinity, being the area averaged values here rice area is dominant in better agreement with observations.
Simulated soil salinity under sesame farming is in good agreement with observed values by Mondal et al. [2001] if salt influx from flooding is included. Mondal et al. [2001] reported the monthly soil salinity values at two villages (Barondanga and Mirzapur) under non-rice crops for the period 1996-1998. These villages are within the Gutudia Union of 60 km 2 total area, from which 93% is protected land. The field for cultivating non-rice crops (sesame) during the winter season relies on groundwater irrigation. Groundwater at both sites were slightly saline with electrical conductivities from 1.60 to 1.99 dS/m at Barodanga and from 1.75 to 2.04 dS/m at Mirzapur. The simulated irrigation water salinity was also saline but with smaller values from 0.49 to 0.88 dS/m. For the simulation, we assumed that sesame is grown from 15th February during 95 days. Simulated soil salinity under sesame farming in the non-protected area for Q8 ensemble ( Figure S3) are in better agreement with observed values than the protected land simulations (Figure 7). Observedand simulated soil salinities were much lower during the wet season (July-November) than the dry season (February-May) at both sites. Observed soil salinities at Barodanga varied from 2.25 to 12.09 dS/m and from 1.52 to 8.74 dS/m at Mirzapur. Simulated soil salinities for the non-protected land varies from 0.82 to 8.54 dS/m for Q8 scenario, from 0.55 to 8.59 dS/m for Q0 scenario and from 1.63 to 11.55 for Q16 scenario. Simulated soil salinities for the protected land were consistently lower than observed values.

Building Confidence on Soil Salinity Model Spatial Variability
The proposed soil salinity model, for the climate ensemble Q0, can accurately reproduce the observed spatial soil salinity variability despite differences between the climate inputs used and the observed temperature and rainfall ( Figure S4). The most comprehensive spatial snap-shot to date on soil salinity in coastal Bangladesh was reported by the Bangladesh Soil Research Development Institute SRDI [2010]. The map shown in Figure Figure 8 does not correspond with the five soil salinity degrees but a combination of them as indicated in the map legend making it not straight forward to compare with the union level simulated results. A comparison of the simulated annual median soil salinity with the area extent affected by different degrees of soil salinity suggest that the Q0 ensemble (i.e., unperturbed physics) is the closest one to the observations (Table 1). Simulated results for Q0 under estimated the area extent of classes 1-4 and over-estimate area extent for class 5. Figure 9 shows the observed annual median soil salinity for years 2001 and 2009 as reported by Dasgupta et al. [2015] and the simulated soil salinity for the different climate ensembles. Again, the simulated results for climate ensemble Q0 are the closest to the discrete observations. The simulated annual median shows the observed gradient of salinity decreasing with distance from the Sundarbans/coast.

Soil Salinity Projections by 2050
Using the proposed simulation model and the Q0 climate scenario, the annual median soil salinity and the third quartile of the daily soil salinity for the dry season (November-April) and wet season(July-September) for year 2050 has been calculated ( Figure 10). The annual median soil salinity captures the central tendency in annual salinity. The third quartile of daily soil salinity data during the dry season captures the median tendency of the upper half of the soil salinity daily values, which permits the assessment of the seasonally intense salinity on important dry season crops (i.e., HYV Boro rice). The third quartile of daily soil salinity data during the wet season permits the assessment of important crops such as paddy-irrigated Aman rice. A visual analysis of the spatial distribution of soil salinity suggest a small change during the wet and dry season but a decrease on the annual median soil salinity at the central area of the study zone that becomes less saline by 2050. The observed (year 2009) and projected (year 2050) annual median and maxima during wet and dry season reported by Dasgupta et al. [2015] are included in Figure 10 for comparison purposes. The annual median soil salinity projections under Q0 scenario are significantly lower (one to two degrees of salinity for more than 60% of the observation stations) than the projections by Dasgupta et al. [2015]. Simulated area extent of soils affected to a certain degree of salinity, based on annual median soil salinity, for climate ensembles Q8 and Q16 are larger than for Q0 (see below) but still systematically smaller than the projected values by Dasgupta et al. (see Figure S5 for Q8 and Q16 results). The simulated soil salinity third quartile is also systematically smaller than the Dasgupta et al.'s projection. These differences might be, to a certain degree, explained by the fact that we are comparing different statistics than reported by Dasgupta et al. (i.e., third quartile of daily values versus maximum of monthly averaged values). We have selected the third quartile because the authors have no access to the time series used by Dasgupta et al. [2015] and because at daily time steps, the maximum soil salinity might represent only a few days that is not comparable to the maximum monthly averaged values reported by Dasgupta.
A comparison of the simulated projections for all three climate ensembles suggest that there is not a clear trend on soil salinity changes by 2050 but is very much controlled by inter-season variability ( Figure 11). Based on annual median soil salinity projections alone, it seems clear that; under a Q0 climate, soil salinity will be better off by 2050 relative to 2009 with an increase of 3,442 km 2 of non-saline soils and a reduction on the area extent of all the saline soil types; under Q8 soil salinity will be on worst off, with a reduction on non-saline soils of 2,500 km 2 and an increase of all types of soil salinity; under Q16 soil salinity will be similar to the one in 2009. Dry season area changes based on the third quartile suggest a different story in which soil salinity changes under Q0 and Q8 are very similar and worst off than 2009 and, again Q16 is very similar to baseline year 2009. Wet season area changes based on the third quartile suggest a similar narrative to the one described for the annual median soil salinity changes but with smaller area changes.

Discussion
We have developed a system dynamic model that integrate for first time the key drivers of soil salinity in low-lying deltas illustrated in Figure 1 at regional scale. Linear emulators allow us to integrate the different hydrological drivers of soil salinity for the large study area without the large computational time penalty. The time required to emulate all the hydrological variables for all 653 unions for 50-year simulation is ∼1 min on a standard 4-Cores 2015 laptop. The accuracy of the emulated daily hydrological time series (river elevation, inundation depth) is well within a factor 10 uncertainty, which is considered adequate for modeling 10.1002/2016EF000530 purposes and tends to be more accurate for the upper range of values which are often more important (e.g., higher river elevations cause tidal flooding and salt ingress). The predictors used for the different emulated variables have been selected following a combination of informed guess and trial and error. The presented predictors are those who best reproduce the simulators and are not over-fitted to accommodate the inherently uncertainty of the simulators. The river area has been subtracted from the total union area to avoid counting rivers as agriculture land. River area for each union is provided in Table S4.
The proposed model structure presented in Section 2 has been able to reproduce the behavior (time series and spatial variability) observed in coastal Bangladesh, building confidence on the validity of the proposed model structure. The approach of using a daily time step and land unit area as proxy for the entire union but divided into protected and non-protected land, seems to capture spatial and temporal variability of soil salinity under different cropping patterns adequately. Dividing the union into protected and non-protected area allows us to obtain a union representative value by area averaging the results. Non-protected areas that  are prone to flooding, tend to have the largest soil salinity values. The daily time step calculation allows capturing not only the seasonal variability but also the rapid soil salinity wash-out during the monsoon season.
We have shown how the model is able to reproduce the seasonal soil salinity increase during the dry season and decrease during the rainy season for integrated rice and shrimp farming and non-rice crops. Using the dominant (e.g., area-wise) cropping pattern for each union, we have shown how the model captures the spatial soil salinity variability for May 2009. The model adequately reproduces the spatial pattern of salinity.
Projections by 2050 of the annual median, and third quartile for the dry and wet season suggest that inter-seasonal climate variability will determine if area of soils affected by different degrees of salinity will be better or worse off than the baseline year 2009. Under the Q0 climate scenario, the area changes based on the annual median soil salinity suggest that an order of 3,500 km 2 of non-saline soils will be recovered. This number is one order of magnitude larger, and of different sign, than the observed increase of 354 km 2 on land affected by various degrees of soil salinity for the 9 years period 2000[SRDI, 2010. Under the Q8 and Q16 climate scenarios, the simulation suggests a net loss of not saline soils in the order of 2,500 km 2 and 270 km 2 , respectively. Projections by 2050 of annual median soil salinity at a number of stations reported by Dasgupta et al. [2015] are systematically larger than the projected soil salinities at union level, and changes relative to baseline year 2009 is of different sign (i.e., soil salinity decreases) for Q0 climate scenario. Clarke et al. [2015] indicates that climatic change up to the year 2100 will cause higher temperatures and hence, increases in crop evapotranspiration, but the inter-annual variability of rainfall and sequences of dry years was found to be the main driver of increased soil salinity. The results presented here reinforce the idea of inter-season variability of rainfall and temperature being one of the main drivers of soil salinity changes.
The proposed method has a number of limitations. The use of emulators to simulate the main hydrological drivers of soil salinity implies that the model is only reliable when the input drivers fall within the range of values used for training the emulators. If any of the simulation data used to train the emulator are biased, the emulator will be also biased. For example, the authors are well aware that the river salinity in the Western region of the study zone is systematically over predicted by the FVCOM simulations. This may be due to missing channels in the delta model such as the Nabaganga river. This channel is not represented in FVCOM, and instead the waters are channeled into the Madhumati (which flows further to the east). There is also the possibility that freshwater sources are missing-as FVCOM is not including evaporation/precipitation. Therefore, soil salinity in this region needs to be interpreted with caution. Our flooding results are realistic, but flooding can change significantly if embankments are breached or elevated. For our flooding simulations, we have assumed that flooding occurs without breaching the embankments (i.e., only overtopping occurs). Groundwater pumping is an input to both the simulated and emulated calculations but are imposed as a time-varying boundary conditions instead of being dynamically linked to the groundwater stock. Due to the large uncertainty on groundwater stock in coastal Bangladesh [Richey et al., 2015], the authors decided not to attempt in this model version to link groundwater depth to changes in calculated irrigation needs. The modulation of the surface drainage rate by the daily river elevation is an over simplification of the complex daily operation of the polders sluice gates. To effect drainage, gates are opened during low tide to allow drainage, and then closed at high tide to prevent ingress of river water (whether saline or not). Gates are also operated in the dry season to allow inflow of river water for irrigation, provided the river water is fresh. The proposed modulated surface drainage rate is not suited to resolve this active sluice gate operation but to capture the inter daily changes on surface drainage due to changes on the river elevation relative to the polder elevation.

Conclusions
Understanding the detailed dynamics of salt movement in the soil is a prerequisite for understanding salinity trends and devising appropriate management strategies to improve land productivity of coastal regions. Numerical models able to resolve soil salinity at regional scale (∼10,000 km 2 ) and the daily variability will allow coastal managers to explore what if case scenarios and proactively manage the land uses for agriculture of low-lying lands such as Deltas. While there are models that resolve the water and salt balance equations, these are either limited to farm land size scale or stream salinity. In this work, we have presented a system dynamic model that resolves the salt and water balance equations at the root zone at regional scale. The proposed approach is demonstrated for the south-west coastal zone of Bangladesh and produces results consistent with the limited observations. Using the proposed approach, we have projected the soil salinity for three different climate ensembles for the year 2050, showing how potential adverse trends varies with the climate scenario. Using these results, we hope to inform farmers, planners and managers on the likely trajectory of salinity development and the impact of future soil salinity on crop production and on farmers' livelihoods in these sensitive coastal ecosystems. (Bangladesh Agricultural Research Institute) and Md. Anwarul Abedin and Abu Zofar Md. Moslehuddin (Bangladesh Agriculture University) for providing observed crop properties and cropping patterns, respectively. This work "Assessing Health, Livelihoods, Ecosystem Services and Poverty Alleviation in Populous Deltas" (NE-J002755-1)' was funded with support from the Ecosystem Services for Poverty Alleviation (ESPA) program. The ESPA program is funded by the Department for International Development (DFID), the Economic and Social Research Council (ESRC) and the Natural Environment Research Council (NERC). The support of the UK British Council INSPIRE R-4 program is appreciated. The authors also acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work. The authors would also like to acknowledge the useful and constructive feedbacks provided by the anonymous reviewer and Dr. Mac Kirby. Following AGU data policy, the authors agree on curate the data used in this study for at least 5 years after publication and provide a transparent process to make the data available to anyone upon request.