Simulation of Salinity Intrusion in Chao Phraya River by a Developed Mathematical Model and MacCormack Scheme with Cubic Spline Interpolation

The purpose of this research was to develop a mathematical model for controlling salinity intrusion in Chao Phraya River, Thailand, using one-dimensional advection-diffusion equation. There have been many research studies that applied a mathematical model called dispersion model to estimating salinity concentration. Our proposed dispersion model is simple, using very few parameters, but can simulate salinity intrusion adequately. We used MacCormack scheme to approximate salinity intrusion and cubic spline interpolation technique to approximate field data including initial salinity concentrations and the concentrations at the left boundary. The MacCormack scheme and Cubic Spline interpolation technique were quite suitable for approximating real data, and the proposed mathematical model was able to predict salinity intrusion adequately.


Introduction
Rivers and canals are vital to the existence and livelihood of people. Chao Phraya River originated from the confluence of the Nan and Ping rivers in Nakhon Sawan Province, Thailand. The river flows southwards for about 370 kilometers through the central plains of the country to Bangkok and then down to the sea at the Gulf of Thailand.
Chao Phraya River is maintained by the Water Transmission and Maintenance Project. A main objective of the project is to fix problems related to drought and salinity intrusion. Chao Phraya Dam releases fresh water down into the river to counter salinity intrusion from the sea [1].
In addition to the lack of fresh water during the dry season, the rise and fall of tides around an estuary is another factor that causes salinity intrusion in the Chao Phraya River. The higher level the tides are, the more salt water is pushed upwards a river. Therefore, the Thai government set a measure to constantly monitor the quality of water sources. The real-time water quality monitoring system is generally updated on an hourly basis (see Figure  1) [2] .
Ten stations along Chao Phraya river monitor the salinity concentration in real-time. If the salinity concentration exceeds the standard, the station will issue a warning to related agencies so that counter measures can be made and implemented (See Figure 2) [3] .
A mathematical model is an inexpensive tool that can help approximate, predict, and manage salinity intrusion related problems. Several researchers have conducted studies related to this issue.   [3] In [4] , two mathematical models-Crank-Nicolson hydrodynamic model and a dispersion model with backward time central space scheme-were used to simulate pollution caused by sewage effluent in the nonuniform water flow in a stream with inconsistency in current velocity. That study demonstrates that the models were applicable to real-world cases. In [5] , a better version of a finite difference scheme was used to solve advection-dispersion-reaction equation (ADRE). In that study, a hydrodynamic model and an advection-dispersion-reaction model were used to simulate pollution from sewage effluent. In [6] and [7] , research has been presented on the increase of salt water.
In [8] , a one-dimensional simulation of water quality 388 Simulation of Salinity Intrusion in Chao Phraya River by a Developed Mathematical Model and MacCormack Scheme with Cubic Spline Interpolation measurement in a stream was done by using a finite volume method to come up with estimated solutions for pollutant concentration equation. The solutions provide logical and rational results that can be applied to water quality measurements. In [9] , an Optimal Homotopy Asymptotic Method or OHAM was used to investigate and simulate dispersion of pollutants in water.
In [10] , the authors focused on salinity intrusion in the area of the Tha Chin estuary, Thailand and showed that their model was suitable for forecasting natural phenomena related to salinity intrusion. In [11] , the researchers came up with a numerical model of salinity's advective transport in the coastal areas, associated with artificial diffusion using upstream differences. This model could be used with small grids of two to four kilometers in the case of tidal problems, while larger grids could be applied to non-tidal problems. In , a tidally and cross-sectionally averaged model were used with Hansen and Rattray equations to simulate the distribution of salinity as well as the vertical exchange flow along the estuary of the Hudson River. Solutions from the model were assessed for the northern San Francisco Bay and the Hudson River over forcing conditions' range.
Some researchers have used Cubic Spline interpolation technique to approximate data related to salinity intrusion and salinity value. In [13] , a connection between cubic splines and a popular compact finite difference formula was detailed. The researcher claimed that the compact approach for nonuniform meshes provide some advantages while the one for uniform meshes provide an additional way that helps treating edge effects. In [14] , two types of cubic spline functions-cubic spline interpolation with 2 C continuity and Piecewise Cubic Hermite Spline (PCHIP) with 1 C continuity-were used for interpolation. In [15] , stream water quality modeling was done by a simple modification of MacCormack and Saulyev schemes. They were used to solve dynamic one-dimensional advectiondispersion-reaction equations (ADRE). The results showed a better prediction accuracy without major loss of efficiency in computation. In [16] , the general form of teeth method was compared with a cubic spline method applied to a mathematical model of teeth of a child in order to obtain a minimum error in showing the form of general teeth while creating a better view of the nature of orthodontic wires. The results of this cubic spline method showed a symmetrical teeth's shape of curve of a normally asymmetrical curve. An example was provided in that study to show how efficient, accurate, and simple the proposed method is.
In this study, we focused on parameter optimization of a mathematical model of Chao Phraya River under the influence of tides that affected salinity values in the river. The salinity concentration generated by the model, MacCormack technique and spline interpolation were compared with real salinity concentration to see whether the model could estimate real data adequately or not.

The Governing Equation
The mathematical model describing transportation and diffusion processes is a one-dimensional advection-diffusion-reaction equation (ADRE).The advection-diffusion equation is [5,10,17 and 18] 2 2 , S x t is the salinity concentration (g/l); U is the water flow velocity in x -directions (m/s); D is the salinity diffusion coefficient (m 2 /s); L is the length of a condimental stream; and T is the stationary time of simulation.
Letting that where s u is the salinity advection rate; w u is the fresh water flow velocity function that is released by the diversion dam; and k is the dilution rate of salinity by fresh water, for all 0 x L < ≤ and 0 .

Initial and Boundary Conditions
The effects of salinity intrusion in this estuarine, which likely will become more severe in the future, are caused by ocean tides. The mathematical model was applied to analyze the salinity intrusion. The initial salinity values were different depending on the time and distance along the river.

The Initial Condition
The initial condition was assumed to be where ( ) f x is a given initial measurement of salinity concentration function along the considered river.

The Boundary Condition
The left boundary condition was assumed to be where ( ) g t is a given present measurement of salinity concentration function at the first monitoring station which is the closest to the estuary. The right boundary condition was assumed to be 0 ( , ) , where 0 k is the rate of change of salinity concentration around the most distant monitoring station.

Numerical Techniques
The salinity concentration was calculated by using MacCormack Scheme, and the values at the left boundary condition and initial values were calculated by using Cubic spline interpolation, as described below.

The First
Step of the MacCormack Scheme The first step of the MacCormack scheme is to approximate S and its derivatives in equation (3) by forward time and forward space scheme (FTFS) as follows: Substituting (7) into (3), we get and 0 1, n N ≤ ≤ − and simplify form (8) by ( ) and ( ) For left boundary, where 1, i = we obtain ( ) ( ) we obtain the MacCormack predictor step formulation ( )

The Second
Step of the MacCormack Scheme The second step of the MacCormack scheme is approximation equation (3) by a backward time and backward space scheme (BTBS) following by: Substituting (17) into, we get and simplify (18)

Cubic Spline Interpolation
Raw salinity concentration data were obtained by field measurement at the first monitoring station which is the closest to the Chao Phraya estuary. Function ( ) f x and ( ) g t were defined by using those data. In this research, cubic spline interpolation was introduced to represent the field data for initial and boundary conditions as described below.
The cubic spline interpolation method, which is a piecewise polynomial approximation, considers a unique cubic equation of each subinterval such that they have necessary conditions as first and second derivatives of a continuous cubic spline [19,20] .
(5). One of sets of boundary conditions is satisfied according to this condition as where T is a tridiagonal matrix of n n × as follows The results of coefficients can always be written in the form of cubic function on [ ]

Numerical Experiment
This study is divided into two main steps. In the first step, hourly-based data collected at 8 monitoring stations along the Chao Phraya River during 05/3-10/2019 were interpolated with cubic spline to the initial condition ( , 0) S x (g/l) and left boundary condition (0, ) S t (g/l), adequately approximating the values at these conditions (see Figure 3 and 4). In the second step, salinity intrusion concentrations at all monitoring stations were numerically simulated and approximated by a simple advection-diffusion equation and MacCormack Scheme. Figure 5 and 6 show the approximated salinity concentrations in a red line together with scattered points of real data from two monitoring stations, one closest to the estuary and another one more distant to the estuary.
The optimum values of s u and w u of our model were found by assigning 6 s u values and 6 w u values (see Figure 7 to 10) to the model and determining the simulated salinity concentration for each combination of these s u and We also used these optimum parameters to run our model at different times and compared the resulting simulated salinity concentrations to the real concentrations at 4, 8, 12, 16, 20 and 24 (hrs.) at one monitoring station and found a good match between those values (See Table 1).

Salinity (g/l)
Comparison between real data and approximated data at a monitoring station Approximation Data

Real Data
The standard level of slinity of raw water to produce water supply should not be over 2.5 g/l

Discussion and Conclusions
In this study, salinity intrusion concentrations along the Chao Phraya River were approximated with a one-dimensional advection-diffusion equation and compared with real concentrations automatically monitored in real-time at several monitoring stations (from which data were affected by tides at the downstream). The values of the parameters ( s u , w u ) of the equation were optimized. The resulting match between the approximated and real values was satisfactory. These optimum parameter values will be used in a future experiment that will test a new, improved model.

Conflicts of Interest
The authors declare that there are no conflicts of interest with any parties whatsoever.