JP2018045661A - Water area environment simulation device, water area environment simulation method, and program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To accurately run a simulation even in the case where a bridge is built in a water area.SOLUTION: A water area environment simulation device comprises: a model storage unit for storing a model of a water flow, set for each grating set in a simulation target water area, including a dominant equation calculating the effects of a bridge pier included in the grating on the water flow with a projection area and volume of the bridge pier as a parameter; and a simulation processing unit for calculating the water flow in the simulation target water area on the basis of the model.SELECTED DRAWING: Figure 1

Description

  The present invention relates to a water area environment simulation device, a water area environment simulation method, and a program.

  Several techniques have been proposed in connection with the simulation of the marine environment. For example, Patent Document 1 describes that meshes are set using a GIS (Geographic Information System) when a simulation of a tidal current flow or the like is performed.

JP 2003-150889 A

  For example, there is a case where a bridge is provided in the water area like Tokyo Bay. In this case, in order to accurately simulate the marine environment, it is necessary to reflect the effect of the pier on the water flow in the simulation.

  The present invention provides a water environment simulation apparatus, a water environment simulation method, and a program capable of performing a simulation with high accuracy even when a bridge is provided in a water area such as the ocean.

  According to the first aspect of the present invention, the water environment simulation device has the effect of the pier on the water flow with the projected area and volume of the pier included in the lattice as parameters for each lattice set in the simulation target water region. A model storage unit that stores a model of a water flow including a governing equation for calculating the water flow, and a simulation processing unit that calculates a water flow in the simulation target water area based on the model.

  The water environment simulation device stores an interstitial distance storage unit that stores an interstitial distance set for each position in the simulation target water area, and lattice point coordinates that calculate coordinates of lattice points based on the interstitial distance And a calculation unit.

  The water area environment simulation device includes a data storage unit that stores data in which a water level, a vertical kinematic viscosity coefficient, and a vertical diffusion coefficient are represented by numerical values having more digits than other numerical values, and the simulation processing unit includes: You may make it calculate the water flow in the water area of the said simulation object using the data which the said data storage part memorize | stores.

  The water environment simulation device stores the number of grids for each row and column of the grid in the horizontal direction, and the row direction within the range of the number of grids stored in the grid number storage unit And a loop processing unit that executes a loop in the column direction, and the simulation processing unit may calculate the water flow in the grid for each grid in the loop executed by the loop processing unit.

  The loop processing unit may execute the calculation of the rows and columns of the grid in the horizontal direction by parallel processing of the same loop.

  The simulation processing unit may calculate water quality and pollution for each grid set in the simulation target water area, and may calculate transparency for each grid point based on the obtained water quality and pollution. .

  According to the second aspect of the present invention, in the water environment simulation method, for each grid set in the simulation target water area, the influence of the pier on the water flow using the projected area and volume of the pier included in the grid as parameters. A water environment simulation apparatus including a model storage unit that stores a model of a water flow including a governing equation for calculating a water flow state calculation step of calculating a water flow in the simulation target water region based on the model.

  According to the third aspect of the present invention, for each grid set in the simulation target water area, the program calculates the influence of the pier on the water flow using the projected area and volume of the pier included in the grid as parameters. A program for causing a computer including a model storage unit that stores a model of a water flow including a governing equation to execute a water flow state calculation step of calculating a water flow in the simulation target water area based on the model.

  According to the present invention, simulation can be performed with high accuracy even when a bridge is provided in the water area.

It is a schematic block diagram which shows the function structure of the water area environment simulation apparatus which concerns on embodiment of this invention. It is explanatory drawing which shows the example of the input-output data with respect to the water area environment simulation apparatus which concerns on this embodiment. It is explanatory drawing which shows schematic structure of the simulation model which the model memory | storage part concerning this embodiment memorize | stores. It is explanatory drawing which shows the example of schematic structure of the elution submodel which concerns on this embodiment. It is explanatory drawing which shows the example of the circulation of the substance relevant to the diatom which the water quality model which concerns on this embodiment simulates. It is explanatory drawing which shows the example of a setting of the variable grating | lattice in this embodiment. It is explanatory drawing which shows the example of the program using the interstitial distance set to the variable of the two-dimensional arrangement | sequence which concerns on this embodiment. It is explanatory drawing which shows the example of a program when not performing an acceleration process. It is explanatory drawing which shows the example of the program at the time of performing a high-speed process in this embodiment. It is explanatory drawing which shows the example of the program which has arrange | positioned the vertical loop inside the horizontal loop in this embodiment. It is explanatory drawing which shows the example of the program which parallelized the loop in this embodiment.

Hereinafter, although embodiment of this invention is described, the following embodiment does not limit the invention concerning a claim. In addition, not all the combinations of features described in the embodiments are essential for the solving means of the invention.
FIG. 1 is a schematic block diagram showing a functional configuration of a water environment simulation apparatus according to an embodiment of the present invention. As shown in the figure, the water environment simulation apparatus 100 includes a display unit 110, an operation input unit 120, a data input / output unit 130, a storage unit 180, and a control unit 190. The storage unit 180 includes a model storage unit 181, an inter-lattice distance storage unit 182, a data storage unit 183, and a lattice number storage unit 184. The control unit 190 includes a simulation processing unit 191. The simulation processing unit 191 includes a lattice point coordinate value calculation unit 192 and a loop processing unit 193.

The water environment simulation apparatus 100 simulates the influence of water currents and living organisms in the ocean, and calculates values indicating the marine environment such as the concentration of a substance that melts in the ocean. The water environment simulation apparatus 100 is configured using a computer.
In the following description, the simulation target is the ocean, but the simulation target of the water environment simulation device 100 is not limited to this. For example, the water environment simulation apparatus 100 may perform simulation of a water area other than the ocean such as a lake.

The display unit 110 includes a display screen such as a liquid crystal panel or an LED (Light Emitting Diode), and displays various images. For example, the display unit 110 displays a simulation result such as the concentration of a substance that melts in the sea.
The operation input unit 120 includes input devices such as a keyboard and a mouse and receives user operations.

The data input / output unit 130 inputs / outputs data to / from the water area environment simulation apparatus 100. For example, the data input / output unit 130 acquires various data for simulation and outputs data indicating the simulation result.
The data input / output unit 130 includes a communication circuit, for example, and communicates with other devices. Alternatively, the data input / output unit may include an interface for a storage device such as a USB (Universal Serial Bus) port, and input / output data to / from the storage device. Alternatively, the data input / output unit 130 may include both a communication circuit and a storage device interface.
Part or all of the simulation data may be input by a user operation from the operation input unit 120. Further, as described above, a part or all of the simulation result may be displayed on the display unit 110 in addition to or instead of the data output from the data input / output unit 130.

The storage unit 180 stores various data. The storage unit 180 may be built in the main body of the water environment simulation device 100. Alternatively, the storage unit 180 may be configured using an external storage device for the main body of the water environment simulation apparatus 100.
The model storage unit 181 stores a model (simulation model) of a simulation performed by the water environment simulation device 100. In particular, the model storage unit 181 includes, for each grid set in the simulation target water area, a water flow including a governing equation that calculates the influence of the pier on the water flow using the projected area and volume of the pier included in the grid as parameters. Remember the model. A lattice here is each cell (Cell) which divided the space (water area) of simulation object into a mesh. Details of the simulation model stored in the model storage unit 181 will be described later.

The interstitial distance storage unit 182 stores the interstitial distance (distance between lattice points) set for each position in the simulation target water area.
The data storage unit 183 stores various data for simulation, such as data measurement values (observation values) in a simulation target water area. In particular, the data storage unit 183 stores data in which the water level, the vertical eddy viscosity coefficient, and the vertical eddy diffusion coefficient are expressed by numerical values having more digits than other numerical values.
The lattice number storage unit 184 stores the number of lattices for each row and column of the lattice in the horizontal direction for the lattice set in the simulation target water area.

  The control unit 190 controls each unit of the water environment simulation device 100 and executes various processes. The control unit 190 is configured by a CPU (Central Processing Unit) included in the water environment simulation device 100 reading out and executing a program from the storage unit 180.

  The simulation processing unit 191 performs a simulation of water flow and water quality using the simulation model stored in the model storage unit 181. In particular, the simulation processing unit 191 simulates the water flow in the simulation target water area based on the simulation model stored in the model storage unit 181. In addition, the simulation processing unit 191 calculates the water flow in the simulation target water area using the data stored in the data storage unit 183. In addition, the simulation processing unit 191 calculates a water flow in the lattice for each lattice in a loop executed by the loop processing unit 193. As a result, the simulation processing unit 191 executes repetitive processing in the simulation. In addition, the simulation processing unit 191 calculates water quality and pollution for each grid set in the simulation target water area, and calculates transparency for each grid point based on the obtained water quality and pollution.

The lattice point coordinate value calculation unit 192 calculates the coordinates of the lattice points based on the interstitial distance stored in the interstitial distance storage unit 182.
The loop processing unit 193 performs loop processing in program execution. In particular, the loop processing unit 193 executes a loop in the row direction and a loop in the column direction within the range of the number of lattices stored in the lattice number storage unit 104 in the simulation performed by the simulation processing unit 191. The loop processing unit 193 executes at least one of the rows and columns of the lattice by parallel processing.

FIG. 2 is an explanatory diagram showing an example of input / output data for the water environment simulation device 100. As shown in the figure, the data input / output unit 130 receives external boundary data, river daily average data, point source daily average data and initial setting data, and the simulation processing unit 191 uses these data. To pre-process.
The external boundary data is data indicating boundary conditions between the water area to be simulated and the external world. The external boundary data includes, for example, the sun-moon synthetic diurnal tide K2, the main lunar half-diurnal tide M2, the main solar half-diurnal tide S2, the main lunar diurnal tide O1, the main solar diurnal tide P1, and the sun-moon synthetic diurnal tide K1. Tide data such as retardation and amplitude of each of the tides, such as Zhou Yuan Long Yen N2.
In addition, the external boundary data is, for example, chemical oxygen demand (COD), total nitrogen (TN), total phosphorus (Total Phosphorus) at the boundary between the water area to be simulated and the external world. TP), including the boundary conditions for salinity and water temperature.

Each river daily average data is data which shows the boundary condition with each river flowing into the simulation target water area as an average value per day. Each river daily average data includes, for example, boundary conditions regarding the flow rate of water, chemical oxygen demand, total nitrogen, and total phosphorus at the boundary between a simulation target water area and a river.
Each point source daily average data is data indicating, for example, an average value per day for a boundary condition in a point source that affects a simulation target water area other than a river, such as drainage from a factory. Each river daily average data includes, for example, boundary conditions regarding the flow rate of water and the like at the boundary between the simulation target water area and the point source, chemical oxygen demand, total nitrogen, and total phosphorus.
The initial setting data is data indicating an initial value in the simulation target water area. The initial setting data includes, for example, data on initial values relating to chemical oxygen demand, total nitrogen, total phosphorus, salinity, and water temperature at predetermined points in the water area to be simulated.

The simulation processing unit 191 performs preprocessing such as conversion processing and distribution processing on these data. For example, the simulation processing unit 191 converts the chemical oxygen demand contained in these data into total organic carbon (TOC) and biochemical oxygen demand (BOD). I do. In addition, the simulation processing unit 191 converts the total nitrogen (concentration of total nitrogen) included in these data into various modes such as suspended organic nitrogen, dissolved organic nitrogen, ammonia nitrogen, nitric acid + nitrite nitrogen ( A distribution process is performed to distribute the density for each mode.
The simulation processing unit 191 stores the data obtained in the preprocessing in the storage unit 180 and uses it for simulation execution (model calculation processing).

The data input / output unit 130 further receives input of topographic data, sea surface wind data, daily average weather data, calculation grid / water depth data, pier resistance data, flow calculation basic conditions, and water quality calculation basic conditions.
The terrain data includes, for example, terrain information such as sea (for example, seafloor terrain), land (especially coastal terrain), rivers (especially estuary position), breakwater information such as breakwater position, landfill location and reclamation time. Including landfill information.

  The sea surface wind data is data indicating the direction and intensity of the wind on the surface (sea surface) of the water area to be simulated. The sea surface wind data includes, for example, information on hourly wind speed (for example, 16 azimuths) and wind speed on the sea surface of the simulation target water area. When the simulation target area is wide, the sea surface wind data may include information for each area obtained by dividing the simulation target water area.

  The daily average weather data is data indicating the weather information in the simulation target water area as an average value per day. The daily average weather data includes, for example, information on wind power at the sea, wind direction (for example, 16 directions), temperature, humidity, cloudiness (for example, 10 minutes ratio), precipitation, and solar radiation. When the simulation target area is wide, the sea surface wind data may include information for each area obtained by dividing the simulation target water area.

The calculation grid / water depth data is data for setting a calculation grid in a simulation target water area. The calculation grid / water depth data includes, for example, information on each grid such as one or both of the four side lengths (interstitial distance) and the center of gravity position of the calculation grid, and the average water depth (information indicating the depth for setting the calculation grid). ).
The interstitial distance storage unit 182 stores the calculation grid / water depth data, and the lattice point coordinate value calculation unit 192 calculates the coordinates of the grid points based on the calculated grid / water depth data stored in the interstitial distance storage unit 182. To do.
As will be described later, the interstitial distance may vary depending on the location.

  The pier resistance data is data for simulating the influence of the structure on the water flow when a structure such as a pier is installed in the simulation target water area. The pier resistance data includes, for example, information indicating the structure position three-dimensionally, such as the center position, shape, and size of the structure, and a resistance coefficient of the structure against running water.

  Here, in the flow calculation, diffusion and water flow are calculated, and in the water quality calculation, the amount of substance carried by the flow is calculated. The value necessary for flow calculation is called flow calculation basic condition. Examples of flow calculation basic conditions include meteorological data, external tide water level data, river flow data, and the like. On the other hand, a value necessary for water quality calculation is called a water quality calculation basic condition. Examples of water quality calculation basic conditions include river inflow material data in addition to flow data, water temperature data, and salinity data obtained by flow calculation.

The simulation processing unit 191 stores these data in the storage unit 180. Then, the simulation processing unit 191 reads the data acquired by the data input / output unit 130 and the data obtained by the preprocessing from the storage unit 180 and executes the simulation.
The simulation processing unit 191 stores the calculation log in the storage unit 180 during execution of the simulation. The calculation log here is information for the user to confirm the state of the water environment simulation device 100 during simulation execution.

The simulation processing unit 191 generates simulation result data such as flow three-dimensional distribution data, ecosystem three-dimensional distribution data, bottom sediment two-dimensional distribution data, suspension state data, and transparency data by executing the simulation. Then, the data input / output unit 130 outputs these simulation result data according to the control of the control unit 190.
The flow three-dimensional distribution data is data indicating the state of seawater in each part of the simulation target water area. For example, the flow three-dimensional distribution data includes data indicating the flow velocity at each point of the calculation grid with respective components (u, v, and w described later) of the three-dimensional coordinates. Further, for example, the flow three-dimensional distribution data includes data indicating salinity concentration and water temperature at each point of the calculation grid.

  Ecosystem three-dimensional distribution data is data indicating the distribution of substances in the sea. For example, ecosystem three-dimensional distribution data can be obtained for each lattice point in the sea, including dissolved oxygen, chemical oxygen demand, total organic carbon, total organic phosphorus, total organic nitrogen, total phosphorus, total nitrogen, suspended organic carbon. , Suspended organic phosphorus, suspended organic nitrogen, dissolved organic carbon, dissolved organic phosphorus, dissolved organic nitrogen, dissolved inorganic phosphorus, ammonia nitrogen, nitric acid + nitrite nitrogen, diatom (carbon), diatom (Phosphorus), diatom (nitrogen), diatom (silica), dinoflagellate (carbon), dinoflagellate (phosphorus), dinoflagellate (nitrogen), persistent organic carbon (land source), persistent Dissolved organic carbon (land source), diatom (population per unit volume), dinoflagellate (population per unit volume), dissolved silica (Available Silica), suspended silica (Particulate Biogenic Silica), zooplankton (Carbon), zooplankton (phosphorus), zooplank Data showing each concentration of nitrogen (nitrogen).

  The bottom sediment two-dimensional distribution data is data indicating the distribution of substances precipitated on the seabed and substances eluted from the seabed into the sea. For example, bottom sediment two-dimensional distribution data includes the amount of phosphate elution, ammonia nitrogen elution, carbon precipitation, nitrogen precipitation, aerobic carbon content, It includes data indicating the nitrogen content of the air layer, the phosphorus content of the anaerobic layer, and the phosphorus content of the anaerobic layer in terms of the amount per unit area. As will be described later, the aerobic layer and the anaerobic layer are layers obtained by dividing the mud layer on the seabed into two.

The suspended state data is data indicating the concentration distribution of suspended solids (SS) in seawater. For example, the simulation processing unit 191 classifies particles that are suspended substances into four levels by size (diameter), and calculates the concentration of each size particle for each calculation grid in seawater. The suspension state data includes data indicating the concentration of each size particle for each of the calculation grids in the seawater.
The transparency data is data indicating the transparency in each part of the water area to be simulated. For example, the transparency data includes data indicating the transparency at each of the lattice points set on the seawater surface of the water area to be simulated.

  Hereinafter, the simulation model will be described. The model storage unit 181 stores this simulation model. And the simulation process part 191 performs a simulation using this simulation model. Hereinafter, the simulation model is also simply referred to as a model.

FIG. 3 is an explanatory diagram showing a schematic configuration of a simulation model stored in the model storage unit 181. As shown in the figure, the simulation model stored in the model storage unit 181 includes a water flow model, a water quality model, an elution sub model, and boundary conditions.
The water flow model is a model for simulating the flow of seawater. Furthermore, the water flow model includes a calculation formula (conservation formula) for water temperature and salinity.

The water quality model is a model for calculating water quality in each part of the simulation target water area. In particular, the water quality model includes a formula for calculating the concentration distribution of organic matter and nutrients in the sea. The nutrient salt here is a salt used for the growth of an organism. Examples of nutrient salts include various aspects of nitrogen, phosphorus, and silica.
The elution submodel is a model for simulating the influence of the interaction between bottom mud and seawater on the concentration distribution of organic matter and nutrients.
The boundary condition indicates, for example, the boundary between the water area to be simulated and the surrounding environment, such as the sea surface, the seabed and the coast, and the influence of the surrounding environment on the water area to be simulated.

In the simulation model stored in the model storage unit 181, a governing equation is provided for each grid that divides the water area to be simulated (that is, for each mesh obtained by mesh division), and the model is represented. The simulation processing unit 191 performs a simulation by performing calculation using a governing equation for each lattice.
In order to obtain the governing equation, three-dimensional orthogonal coordinates of the x-axis, y-axis, and z-axis are set in the water area to be simulated. The x-axis has an east as an increasing direction, the y-axis has an increasing direction as the north, and the z-axis has an upward direction as an increasing direction. Further, a three-dimensional mesh is set in the water area to be simulated, and each lattice point of the mesh is indicated by three-dimensional coordinates of x, y, and z. The free water surface is expressed as η (x, y, t) using the function η. Here, t indicates time. In addition, the sea floor is expressed as -H (x, y) using the function H.
The simulation processing unit 191 describes a model with simultaneous equations (dominant equations) for each lattice point, and solves the simultaneous equations for each lattice point at each sampling time. The governing equation uses Hydrostatic assumptions and Boussinesq approximation.
Hereinafter, the scalar product is represented by either “•” or “×”. A scalar product is also indicated by parallelizing the objects to be multiplied without expressing an operator.

<Water flow model>
The water flow model includes a continuity equation and a momentum equation for calculating the flow of seawater, a conservation equation for simulating turbulent flow, and a conservation equation for simulating the flow of water temperature and salinity.
(Continuous equation)
The continuity equation is shown as in equation (1).

Here, u, v, and w are ensemble average flow velocities corresponding to the x, y, and z directions, respectively. Hereinafter, the ensemble average flow velocity is also simply referred to as a flow velocity.
(Momentum equation)
The equation of motion relating to the flow velocity u is shown as in equation (2).

As described above, t indicates time. Moreover, ρ ₀ indicates the basic density of seawater. The basic density here is a density at a reference pressure (for example, 1 atm). P represents a pressure (water pressure). K _M represents the vertical eddy viscosity coefficient. f represents a Coriolis parameter, and f _u , f _v , and f _w represent a u component, a v component, and a w component of the Coriolis parameter f, respectively.
Further, F _u is a term indicating a motion generated by a small-scale process that cannot be directly expressed by the grid (subgrid scale) of the model, and is expressed as Equation (3).

Here, A _M is a constant representing the horizontal vortex viscosity. Method will be described later determination of the value of A _M.
In equation (3), the motion generated by the small-scale process is expressed in a representation similar to turbulent diffusion. The same applies to Formula (5), Formula (8), Formula (10), and Formula (12) described later.
Further, the equation of motion relating to the flow velocity v is shown as in equation (4).

Here, F _v is a term indicating the movement caused by the process of direct not express small-scale by the model of the grid, is represented by equation (5).

  Further, the equation of motion relating to the flow velocity w is expressed as in equation (6).

Here, ρ represents the density of seawater. The density ρ of seawater will be described later with reference to equation (16). g represents gravitational acceleration.
(Water temperature and salinity preservation formula)
The conservation equation for the water temperature θ is expressed as in Equation (7).

Here, K _H is a constant that indicates the vertical vortex diffusion coefficient. _Fθ is a term indicating a motion caused by a small-scale process that cannot be directly expressed by the model grid, and is expressed as in Equation (8).

Here, A _H is a constant indicating the horizontal eddy diffusion coefficient. Method will be described later determination of the value of the A _H.
Moreover, the preservation | save formula of the salt concentration S is shown like Formula (9).

Here, F _S is a term indicating a motion generated by a small-scale process that cannot be directly expressed by the grid of the model, and is expressed as in Expression (10).

  Further, the diffusion equation relating to the water quality item M is expressed as in Expression (11). Specifically, M represents the concentration of a substance contained in seawater.

Here, F _M is a term indicating the movement caused by the process of the small-scale that can not be directly expressed by the model of the grid, is represented by equation (12).

(Conservation formula for turbulent energy and turbulent scale)
If there are plant communities etc. in the flow field, it can affect running water and substance mixing. In these cases as well, the individual elements, shapes, locations, etc. in the plant community are irregular and small in scale compared to the spatial scale that needs to be handled as a whole, so they are not expressed individually. This is expressed by setting parameters such as a resistance coefficient. Individual elements locally change the flow, generate turbulence, and dissipate its energy (turbulent energy). The influence of this turbulent energy on material mixing is shown by the following two equations.
The first equation is shown as equation (13).

Here, q ^denotes a turbulent energy at ^q 2/2.
The vector v indicates iu + jv. However, i shows a real axis and j shows an imaginary axis. Therefore, the real component of vector v is u and the imaginary component of vector v is v. u and v are the ensemble average flow velocities described above. The vector v is indicated by an arrow (“→”) above “v”. In Expression (13), “·” indicates an inner product of vectors.
K _q represents the vertical diffusion coefficient (the vertical diffusion coefficient) of turbulent energy.
B ₁ represents _one of the turbulent closure constants. Specifically, B ₁ is a constant obtained empirically in a state where generation and attenuation of turbulent energy are balanced, and for example, B ₁ = 16.6.
Moreover, l shows a turbulent space scale (turbulent length scale). The value of l is determined experimentally or empirically depending on the flow conditions.
F _q indicates the horizontal diffusion of turbulent energy.
Further, the second equation is expressed as shown in Equation (14).

Here, ql indicates a turbulent flow scale (scale of turbulent flow length).
E ₁ and E ₂ are both parameters for matching the model to the observed value. The values of E ₁ and E ₂ are determined experimentally based on observed flow rates. For example, E ₁ = 1.8 and E ₂ = 1.33.
Κ represents the Kalman constant. F _ql indicates turbulent scale horizontal diffusion.
Further, L is a scale representing the distance from the sea surface to the sea bottom, and is expressed as in Expression (15).

Here, η is a function η (x, y, t) indicating the above-described free water surface. In other words, η in equation (15) represents the z coordinate value (height from the average water level (z = 0)) of the sea surface. Further, h in Expression (15) indicates the z coordinate value of the seabed.
Equation (13) shows that turbulent energy is generated by the generation of turbulent flow due to the frictional resistance of the average flow, attenuation of turbulent energy due to the conversion from turbulent energy to thermal energy, and advection and diffusion of the average flow of turbulent energy. It shows that the temporal and spatial distribution is determined. Equation (13) is a so-called turbulent energy prediction equation.
Equation (14) is a governing equation for predicting the turbulent space scale. The turbulent macro scale l represents the largest scale of the turbulent flow, and is an important quantity for determining the attenuation of the turbulent energy and the diffusion coefficient.
The density ρ of the seawater is expressed as in equation (16) by the equation of state.

That is, the density ρ of seawater is represented by a function having the water temperature θ and the salinity concentration S as parameters.
Further, the vertical vortex viscosity coefficient _{K M} and the vertical vortex diffusion coefficient _{K H} is sequentially calculated in accordance with Equation (17).

Here, l is the turbulent macro scale as described above. q indicates the turbulent energy at q ^2/2 as described above. S _M and _SH are stable functions. S _M and _SH are expressed as a function of, for example, the Blacks Richardson number R _f . l and q are determined experimentally or empirically.
In inner bays and coastal areas, structures such as piers may exist in the flow field. In this case, a turbulent flow is generated in the flow field by the structure. Therefore, instead of the equations (13) and (14), a turbulent energy equation and a turbulent scale equation to which terms F _Dq and F _Dql indicating the influence of resistance by the pier are added are used.
The turbulent energy equation to which the term F _Dq is added is shown as equation (18).

The term F _Dq is shown as in equation (19).

Here, C _fq is an experimental constant (a constant determined by experiment). Thus, by adding the term _FDq indicating the influence of the resistance by the pier to the equation (18) indicating the turbulent energy equation, the influence of the structure on the flow field can be expressed more accurately.
Further, F _Dx is expressed as in Expression (20).

Here, C _d denotes a resistance coefficient of the flowing water by piles.
Moreover, A shows the projection area of the pile in 1 lattice. The projected area of a pile within one grid is the product of the diameter and total length of the piles that fall within one grid. For example, in the case of a jacket structure composed of main piles and slant piles having different thicknesses, the projected area A of the piles in one lattice is obtained by Expression (21).

Δx represents a lattice interval in the x direction (x axis direction). D indicates the layer thickness (vertical lattice spacing). Here, in the water environment simulation apparatus 100, a multilayer level model (a model in which a plurality of grids are set in the vertical direction) is used, and the simulation calculation unit 191 calculates a governing equation for each grid constituting each layer. Go and calculate the simulation results. Therefore, the simulation calculation unit 191 obtains the projected area and volume of the pile for each projected area and volume of the pile in one grid (that is, for each grid of each layer).
_CM represents an additional mass coefficient. Here, when a pillar is erected in the sea area, there is a pressure gradient associated with the shape resistance and the acceleration motion of water. A force proportional to the volume obtained by multiplying the pressure gradient accompanying the acceleration motion of water by the volume acts on this column. This factor is referred to as added mass coefficient C _M. For example, in the case of a cylinder, it is possible to use 2.0 as the value of the load weight counting C _M.
V represents the volume of one lattice.
Also, V _d represents the volume of pile in the first grating. The volume of piles in one grid is the product of the cross-sectional area and the total length of the piles that enter one grid. For example, in the case of a jacket structure composed of main piles and slant piles having different thicknesses, the volume V _d of the piles in one lattice is obtained by Expression (22).

Further, _FDy in the equation (19) is expressed as in the equation (23).

_Further , the turbulent flow scale equation to which the term F _Dql is added is shown as the equation (24).

The term F _Dql is shown as in equation (25).

Thus, by adding the term F _Dql to the equation (24) indicating the turbulent scale equation, the influence of the structure on the flow field can be expressed more accurately.
The value of the horizontal vortex viscosity _{A M} can be calculated according to equation (26).

Here, Δy indicates the lattice spacing in the y direction (y-axis direction).
The value of the horizontal vortex diffusion coefficient _{A H} is determined according to equation (27).

The values of constants C _M and C _H are set according to the size of the lattice. For example, C _M = 0.01 and C _H = 0.01 are set for a grid of 1 km (kilometers) × 1 km.

<Boundary conditions>
(Formula of boundary condition on free water surface)
The boundary conditions on the free water surface (sea surface) are set as in the following first to third formulas.
The first equation of the boundary condition on the free water surface is shown as equation (28).

Here, (τ _0x , τ _0y ) indicates the stress due to wind on the water surface. The values of (τ _0x , τ _0y ) are set in advance based on, for example, weather data and theoretical values.
The second equation of the boundary condition on the free water surface is expressed as equation (29).

Here, H ′ represents a pure heat flux on the water surface. S ′ represents the salinity flux on the water surface. The flux here is a flux, that is, an amount (amount of movement) flowing per unit time and unit area.
The pure heat flux H ′ on the water surface is expressed as in equation (30).

Here, phi _S denotes solar radiation amount reaching the earth's surface from the sun (solar radiation). φ _a represents the amount of atmospheric radiation that reaches the ground surface after sunlight is absorbed and scattered by air molecules, aerosols, cloud particles, and the like in the atmosphere. φ _br indicates the amount of infrared radiation that the ground surface emits as a black body. The values of φ _S , φ _a and φ _br are set in advance based on, for example, weather data (meteorological observation data or weather estimated values) and theoretical values.
Φ _Sr indicates the amount of solar radiation (irradiation amount) reflected by the ground surface, and is expressed as in equation (31).

Here, φ _Sc indicates the amount of solar radiation in fine weather. C represents the cloud amount. The values of φ _Sc and C are preset based on, for example, weather data (meteorological observation data or weather estimated values) and theoretical values.
Φ _ar represents the amount of atmospheric radiation reflected from the ground surface, and is expressed as in equation (32).

Here, e indicates the water vapor pressure. T _a indicates the temperature. The value of e and T _a is, for example, preset based on meteorological data (meteorological data or weather estimates) and the theoretical value.
Φ _e represents the latent heat transport amount accompanying evapotranspiration, and is represented by the equation (33).

Here, W _z indicates the wind speed at a height z on the water surface. e _s denotes a saturated vapor pressure for the temperature on the water surface. e _z represents the water vapor pressure at sea level on a height z. T _S indicates the water temperature of the water surface.
Also, phi _c represents the convective heat transfer amount, is represented by equation (34).

Here, T _z indicates the air temperature at a height of Z meters (m) above the water surface. For example, values at a height of about 1 to 2 m above the water surface are used as the temperature T _z and the wind speed W _z . The values of W _z , e _s , e _z , T _S and T _z are preset based on, for example, weather data (meteorological observation data or weather estimated values) and theoretical values.
The third equation of the boundary condition on the free water surface is expressed as equation (35).

  Here, the salinity flux S ′ on the water surface is represented by the equation (36).

Here, [E′−P ′] indicates a pure freshwater flux on the water surface due to evaporation / rainfall. S (0) indicates the salinity concentration on the water surface. The value of [E′−P ′] is set in advance based on, for example, weather data (meteorological observation data or weather estimated value) and a theoretical value. The value of S (0) is determined in advance based on, for example, an actual measurement value.
(Boundary condition equation at the sea floor)
In addition, boundary conditions are set as follows for stress, pure heat flux, and salinity flux at the seabed. At the sea floor, z = −H. In addition, the vertical gradient of θ and S is zero on the sea floor, and there is no heat transfer or salt transfer through the boundary.
The first equation of the boundary condition at the sea floor is shown as equation (37).

Here, (τ _bx , τ _by ) indicates the stress at the seabed. This (τ _bx , τ _by ) is given by a logarithmic approximation law as shown in equation (38).

Here, _{C D} represents the drag coefficient. The value of the drag coefficient _{C D} is 0.0025, or, using the larger value one of the value according to formula (39).

Here, H (x, y) represents the seabed topography (the value of z on the seabed). H (x, y) is also simply expressed as H. Z _b indicates a calculation grid closest to the seabed. u _b represents the flow velocity in the x direction in the calculation grid z _b . v _b represents the flow velocity in the y direction on the calculation grid z _b .
Further, ln represents a natural logarithm. The value of z ₀ depends on the local ancestry of the seabed. For example, the value of z ₀ is 1 centimeter (cm).
Further, the second equation of the boundary condition at the sea floor is expressed as equation (40).

  The third equation of the boundary condition at the sea floor is expressed as equation (41).

Here, Es indicates the amount of sedimentation of the substance particles.
In addition to the free water surface and sea bottom described above, rivers, point sources for draining water, and the like can be modeled under boundary conditions. The boundary condition may be indicated by a numerical value such as a constant or a function value, may be indicated by an equation, or may be indicated by a combination thereof.

<Coordinate transformation>
The Cartesian coordinates (x, y, z, t) are not effective for expressing rapid changes in water depth along the coast as in the inner bay. Therefore, new coordinates (x ^* , y ^* , z ^* , t ^* ) are introduced in the vertical direction for each equation set as described above, and coordinate conversion of all equations is performed. Specifically, the conversion is performed as shown in Expression (42).

  Here, D = H + η. Therefore, σ = 0 when z = η (that is, the water surface). In addition, σ = −1 at z = −H (that is, the sea floor).

<Elution submodel>
The ecosystem model incorporates the interaction between bottom mud and seawater, which is an important process in closed and shallow waters, as an elution submodel. For example, as elements of the elution submodel, suspended organic carbon, suspended organic nitrogen, NH ₄ , NO ₃ , suspended organic phosphorus, PO ₄ , silicon (silica), hydrogen sulfide (sea water), and dissolved oxygen It can be employed, but is not limited to this.
An existing model such as a model developed by EPA (United States Environmental Protection Agency) can be used as an elution submodel used in the water environment simulation apparatus 100.

  In the model, the bottom mud is divided into an aerobic layer and an anaerobic layer. The upper layer becomes either an aerobic layer or an anaerobic layer depending on the dissolved oxygen concentration of the directly above water. On the other hand, the lower layer is always an anaerobic layer. The organic matter settles in the bottom mud from the water, and the mineralization (decomposition) of the organic matter in the bottom mud is shown at three speeds (easily degradable, hardly degradable, and nondegradable). This decomposition shows that organic matter returns to nutrients and oxygen demand (SOD) occurs in the process. In addition, the reverted nutrients are eluted in seawater.

FIG. 4 is an explanatory diagram illustrating an example of a schematic configuration of an elution submodel. As shown in FIG. 4, the elution submodel shows the movement of organic matter, nutrient salts, and dissolved oxygen between the bottom mud and seawater. Further, as described above, the bottom mud is divided into an aerobic layer and an anaerobic layer. The elution submodel shows the transfer of organic matter and nutrients between the aerobic and anaerobic layers. In addition, the elution submodel shows mineralization (decomposition) from organic matter to nutrient salts and deposition as a deposit.
Hereinafter, the elution submodel will be described with reference to mathematical expressions.

(organic matter)
The decomposition of organic substances is classified into three types of decomposition rate classes. G ₁ is a class that halves in 20 days, G ₂ is a class that halves in 1 year, and G ₃ is a class that does not undergo decomposition until it settles into a deposited layer.
The increase / decrease in the particulate organic carbon in the aerobic layer is represented by the formula (43).

Here, “(T-20)” indicates that the standard water temperature in the bottom mud elution rate test is 20 ° C.
The increase / decrease in the particulate organic nitrogen in the aerobic layer is represented by the formula (44).

  The increase / decrease in the particulate organic phosphorus in the aerobic layer is shown by the formula (45).

Here, H ₁ indicates the thickness of the aerobic layer. The value of an H ₁ is for example 0.1 centimeters (cm).
Further, i takes one of 1, 2, and 3. i values ₁ , ₂ , and ₃ are associated with classes G ₁ , G ₂ , and G ₃ , respectively. POC ₁ , PON ₁ , POP ₁ , _{and 1} respectively represent the particulate organic carbon concentration, nitrogen concentration, and phosphorus concentration (unit: mg / l) in the aerobic layer. POC _{, 0} , PON _{, 0} , POP _{, 0} indicate the particulate organic carbon concentration, nitrogen concentration and phosphorus concentration in seawater, respectively. POC 1 _{, 2} , PON ₂ , POP _{2, 2} indicate the particulate organic carbon concentration, particulate organic nitrogen concentration, and particulate organic phosphorus concentration in the anaerobic layer, respectively. Note that POC, PON, and POP represent the particulate organic carbon concentration, the particulate organic nitrogen concentration, and the particulate organic phosphorus concentration, respectively.
_{Further, f POC, i, f PON} , i, f POP, i denotes each particle state organic carbon concentration, particle state organic nitrogen concentration, the distribution coefficient to the class _{G i} of particle states organic phosphorus concentration. By multiplying the density by the distribution coefficient, the density for each class can be obtained. The distribution coefficient is given as a constant in advance, for example. Since the sedimentation speed varies depending on the size of the particles, the calculation can be performed with higher accuracy by dividing each class.

For example, f _{POC, 1} = 0.65, f _{POC, 2} = 0.20, f _{POC, 3} = 0.15.
For example, f _{PON, 1} = 0.65, _{fPON, 2} = 0.28, _{fPON, 3} = 0.07.
Further, for example, f _{POP, 1} = 0.65, f _{POP, 2} = 0.255, f _{POP, 3} = 0.095.

K _{POC, i} , k _{PON, i} , k _{POP, i} indicate the decomposition rates of POC 1 _{, i 2} , PON 1 _{, i} POP 1 _{, i} , respectively. The values of k _{POC, 1} , k _{PON, 1} and k _{POP, 1} are set to 0.035 (unit: day ^−l ), for example.
Further, the decomposition rate of POC ₂ , PON ₂ , POP _{2, 2} is set to 0.0018 (unit: day ⁻¹ ), for example.
Further, the decomposition rate of POC ₃ , PON ₃ , POP ₃ is set to 0.0 (unit: day ^−l ), for example.

In addition, θ _{POC, i} , θ _{PON, i} , θ _{POP, i} indicate decomposition temperature coefficients, respectively.
For example, the values of θ _{POC, 1} , θ _{PON, 1} , θ _{POP, 1} are all set to 1.100.
For example, the values of θ _{POC, 2} , θ _{PON, 2} and θ _{POP, 2} are all set to 1.150.

For example _, the values of θ _{POC, 3} , θ _{PON, 3} , θ _{POP, 3} are all set to 0 (none).
W ₂ represents the deposition rate (unit: m / day (meter / day)). W _net indicates the sedimentation rate from seawater to bottom mud (unit: m / day). W ₁₂ represents a mixed mass transfer coefficient. Deposition rate W ₂ and sedimentation rate W _{net Non} from seawater to sediment is preset, for example, constant. The calculation of the mixing mass transfer coefficient W _12, will be described later with reference to equation (58).

The first term on the right side of the formula (43) indicates a decrease in the particulate organic carbon concentration due to decomposition. The second term shows a decrease in the particulate organic carbon concentration due to sedimentation from the aerobic layer to the anaerobic layer. The third term shows an increase in the particulate organic carbon concentration due to sedimentation from seawater to the aerobic layer. The fourth term shows an increase in the particulate organic carbon concentration by mixing the aerobic layer and the anaerobic layer.
The terms on the right side of the equations (44) and (45) are also the same as in the equation (43) by replacing the particulate organic carbon concentration with the particulate organic nitrogen concentration and the particulate organic phosphorus concentration, respectively.

  Moreover, increase / decrease of the particulate organic carbon in an anaerobic layer is shown like Formula (46).

  Moreover, increase / decrease of the particulate organic nitrogen in an anaerobic layer is shown like Formula (47).

  In addition, the increase or decrease in particulate organic phosphorus in the anaerobic layer is represented by the formula (48).

Here, H ₂ represents the thickness of the anaerobic layer. The value of H ₂ is, for example, 10 centimeters (cm).
The first term on the right side of the formula (46) indicates a decrease in the particulate organic carbon concentration due to decomposition. The second term shows a decrease in particulate organic carbon concentration due to deposition as a deposit. Since organic matter and the like deposited as a buried material is considered not to diffuse into the aerobic layer or the sea, it is excluded from the concentration in the anaerobic layer. The third term of the formula (46) indicates an increase in the particulate organic carbon concentration due to the sedimentation from the aerobic layer to the anaerobic layer. The fourth term shows a decrease in the particulate organic carbon concentration due to an increase in the aerobic layer and the anaerobic layer.
The terms on the right side of the equations (47) and (48) are the same as those in the equation (46) by replacing the particulate organic carbon concentration with the particulate organic nitrogen concentration and the particulate organic phosphorus concentration, respectively.
Further, the first equation of the conservation equation is expressed as equation (49).

  Further, the second equation of the conservation equation is expressed as equation (50).

Here, C _T0 , C _T1 , and C _T2 indicate the concentrations of organic substances. The subscript values “0”, “1”, and “2” indicate seawater, aerobic layer, and anaerobic layer, respectively. K _L12 indicates a diffusion coefficient between the aerobic layer and the anaerobic layer. f _d0 indicates the dissolved state fraction (dissolved fraction coefficient) of the sea area.
κ ₁ ² will be described later with reference to equation (72). κ ₂ indicates the reaction rate (unit: m / day) of the aerobic layer. kappa ₂ values are set to, for example, 0.0.
Further, in Formula (49) and (50), calculated divided into a particle state fraction _{f p} Dissolved fraction _{f d.} Dissolved fraction f _d1 aerobic layer is shown as equation (51).

Particle state fraction f _p1 aerobic layer is shown as equation (52).

The dissolved fraction f _d2 of the anaerobic layer is represented by the formula (53).

The particle state fraction _fp2 of the anaerobic layer is represented by the formula (54).

Here, m ₁ and m ₂ represent the particle concentration in the aerobic layer and the anaerobic layer, respectively. π ₁ and π ₂ represent distribution coefficients in the aerobic layer and the anaerobic layer, respectively. The values of m ₁ , m ₂ , π ₁ and π ₂ are set in advance as constants, for example.
Moreover, the mass transfer coefficient s between seawater and bottom mud is shown as in equation (55).

Here, SOD shows the oxygen consumption rate by bottom mud. [O2 (O)] represents oxygen in water.
J _T1 is expressed as shown in Equation (56).

J _T2 is expressed as in Expression (57).

Here, _{kPOM, i} indicates the decomposition rate of suspended organic matter. The values “1” and “2” of i are associated with classes G ₁ and G ₂ , respectively. θ _{POM, i} ^(T-20) represents the temperature coefficient of decomposition. 20 ° C. is set as the standard temperature. POM _i1 and POM _i2 indicate the concentrations of suspended organic substances in the aerobic layer and the anaerobic layer, respectively.
The mixing mass transfer coefficient _{W 12} is shown as equation (58).

Here, _Dp represents a particle mixing diffusion coefficient (unit: m ² / day). The value of the particle mixture diffusion coefficient _{D p} is set to, for example, 1.2E-04 (floating-point representation).
Further, theta _Dp denotes the temperature coefficient of the _{D p.} The value of the temperature coefficient θ _Dp is set to 1.117, for example.
G _{POC, R} indicates a POC reference concentration (unit: mgC / g) at a water temperature of 20 ° C. in the aerobic layer. The value of G _{POC, R} is set to 0.1, for example.
_{KM and Dp} represent the particle mixing half-saturation constant (unit: mg / L) with respect to oxygen. The values of _{KM and Dp} are set to 4.0, for example.
The mixing of dissolved substances between the aerobic layer and the anaerobic layer becomes larger than the molecular diffusion due to the mixing accompanying the activity of bentos. This effect is incorporated by equation (59).

Here, D _d represents a pore water diffusion coefficient (unit: m ² / day). The value of D _d is set to, for example, 1.0E-3 (floating point representation).
Θ _Dd represents the temperature coefficient of Dd. The value of θ _Dd is set to 1.08, for example.

(ammonia)
The elution of ammonia in the aerobic layer is represented by the formula (60).

Here, NH ₄ (0), NH ₄ (1), and NH ₄ (2) indicate ammonia concentrations in seawater, an aerobic layer, and an anaerobic layer, respectively. That is, 0, 1, and 2 in parentheses indicate seawater, aerobic layer, and anaerobic layer, respectively. As described above, f _d0 represents the dissolved state fraction coefficient of the sea area.
The elution of ammonia in the anaerobic layer is represented by the formula (61).

Here, J _N1 is expressed as shown in Expression (62).

J _N2 is expressed as in Equation (63).

Regarding ammonia, Formula (64) is used instead of Formula (51) as a formula representing f _d1 .

Here, NH ₄ indicates the concentration of ammonia. Equation (52) is used for f _p1 . For ammonia, the value of f _d1 in equation (52) is based on equation (64).
Regarding ammonia, Formula (65) is used instead of Formula (53) as a formula representing f _d2 .

using equation (54) for f _p2. For ammonia, the value of f _d2 in equation (54) is based on equation (65).
Further, the nitrification action received by ammonia is shown by the equation (66).

Here, _κNH4,1 indicates a nitrification reaction rate (unit: m / day). The value of κ _NH4,1 is set to 0.131, for example.
Θ _NH4 represents the temperature coefficient of the nitrification reaction. The value of θ _NH4 is set to 1.123, for example.
KM _{, NH4} represents the nitrification reaction half-saturation constant (unit: mgN / m ³ ) of ammonia. The values of _{KM and NH4} are set to 728, for example.
Further, θ _{KM, NH 4} indicates the temperature coefficient of the nitrification reaction half-saturation constant. The values of _{θKM, NH4} are set to 1.125, for example.
KM _{, NH4, and O2} represent a nitrification reaction half-saturation constant (unit: mgO ₂ / L) for oxygen. The values of _{KM, NH4, and O2} are set to 0.37, for example.

Further, π ₁ NH ₄ represents a fraction coefficient (unit: L / kg) of the aerobic layer. The value of π ₁ NH ₄ is set to 1.0, for example.
Further, π ₂ NH ₄ represents a fraction coefficient (unit: L / kg) of the anaerobic layer. The value of π ₂ NH ₄ is set to 1.0, for example.
J _N1 represents a decomposition nitrogen flux (unit: Mg · N / m ² / day) of the aerobic layer. The value of J _N1 is set to 0.0, for example. J _N2 represents a decomposition nitrogen flux (unit: Mg · N / m ² / day) of the anaerobic layer. The J _N2 with _{J N} referred to. Parameters are set according to conditions such as the sea area and the calculation gauge.
In order to clarify the notation, variables may be enclosed in []. For example, in formula (66), “NH ₄ (1)” is enclosed in [] for the sake of clarity.

(nitric acid)
The elution of nitric acid in the aerobic layer is represented by the formula (67).

Here, NO ₃ (0), NO ₃ (1), and NO ₃ (2) indicate the concentrations of nitric acid in seawater, an aerobic layer, and an anaerobic layer, respectively.

  The elution of nitric acid in the anaerobic layer is represented by the formula (68).

Here, J _NO31 is expressed as shown in Equation (69).

Here, J [NH ₄ ] represents an elution flux into seawater.
With respect to the nitrate, kappa ₁ ² is shown as equation (70).

Further, with respect to nitric acid, κ ₂ is represented by the formula (71).

Here, _κNO3,1 indicates the denitrification reaction rate (unit: m / day) of the aerobic layer. The value of _κNO3 , 1 is set to 0.10, for example.
Further, _κNO3,2 represents a denitrification reaction rate (unit: m / day) of the anaerobic layer. The value of _κNO3 , 2 is set to 0.25, for example.
Θ _NO3 represents the temperature coefficient of denitrification. The value of θ _NO3 is set to 1.08, for example.
Regarding nitric acid, formula (72) is used instead of formula (51) as a formula representing f _d1 .

Equation (52) is used for f _p1 . For nitric acid, the value of f _d1 in equation (52) is based on equation (72).
Regarding nitric acid, Formula (73) is used instead of Formula (53) as a formula representing f _d2 .

using equation (54) for f _p2. For nitric acid, the value of f _d2 in equation (54) is based on equation (73).

(Sulfide)
The elution of sulfide in the aerobic layer is represented by the equation (74).

Here, S (0), S (1), and S (2) indicate sulfur concentrations in seawater, an aerobic layer, and an anaerobic layer, respectively.
The elution of sulfide in the anaerobic layer is represented by the formula (75).

Here, regarding sulfur, Formula (76) is used instead of Formula (51) as a formula representing f _d1 .

Equation (52) is used for f _p1 . For sulfur, the value of f _d1 in equation (52) is based on equation (76).
As for sulfur, Formula (77) is used instead of Formula (53) as a formula representing f _d2 .

using equation (54) for f _p2. For sulfur, the value of f _d2 in equation (54) is based on equation (77).
Here, π _1S and π _2S indicate distribution coefficients. The values of π _1S and π _2S are set in advance as constants, for example.
Here, J _s2 is expressed as in Expression (78).

Α _{O2 and NO3} indicate oxygen consumption (unit: g O ₂ / m ³ ) due to denitrification. The value of αO2 _{, NO3} is set to 2.8571, for example. α _{O2, C} indicates the amount of organic carbon consumed by denitrification. JC indicates the carbon flux to the sea.
In addition, with respect to sulfur, κ ₁ ² is represented by the formula (79).

Here, κ _{H2S, d1} indicates the oxidation reaction rate (unit: m / day) of the dissolved sulfide. The value of κ _{H2S, d1} is set to 0.20, for example.
Moreover, (kappa) _{H2S, p1} shows the oxidation reaction rate (unit: m / day) of a particulate sulfide. The value of κ _{H2S, p1} is set to 0.40, for example.
Θ _H2S represents a temperature coefficient of sulfide oxidation. The value of theta _{H2 S} is set to, for example, 1.08.

_{KM, H2S, and O2} represent sulfide oxidation constants (unit: mg O2 / L). The values of _{KM, H2S, and O2} are set to 4.0, for example.
Further, π _1S represents the aerobic fractionation coefficient (unit: L / kg). The value of π _1S is set to 100, for example.
Further, π _2S represents an anaerobic fractionation coefficient (unit: L / kg). The value of π _2S is set to 100, for example.

(Dissolved oxygen)
The oxygen consumption amount SOD (unit: g O ₂ / m ³ ) is the oxygen demand amount (CSOD, unit: g O ₂ / g N) due to methane oxidation and the oxygen consumption amount associated with nitrification as shown in the equation (80). It is shown as the sum of (NSOD).

  Moreover, oxygen demand CSOD accompanying organic carbon decomposition in bottom mud is shown like Formula (81).

  Moreover, oxygen consumption NSOD accompanying the nitrification action in bottom mud is shown like a formula (82).

Here, [ΣH ₂ S (1)] indicates the total sulfide concentration (unit: g O2 / m3) of the aerobic layer.
Further, α _{O2 and NH4} indicate oxygen consumption (unit: g O2 / gN) due to nitrification. The value of αO2 _{, NH4} is set to 4.5714, for example. θ _NH4 ^(T-20) represents a temperature constant of nitrification at a water temperature of 20 ° C. As the value of θ _NH4 ^(T-20) , for example, a constant in the range of 1.076 to 1.127 can be used. k _NH4,1 ² shows the values of digestibility constants in aerobic layer. As k _NH4,1 ² values, for example, it can be used 0.894. k _{M, NH4} represents a half-saturation constant of nitrification. As the values of k _{M and NH 4} , for example, constants in the range of 0.63 to 1.19 can be used. θ _{M, NH4} ^{(T = 20)} indicates a temperature constant at a water temperature of 20 ° C. K _{O2 and NH4} indicate a half-saturation constant of oxygen concentration in the anaerobic layer. As the value of K _{O2 and NH4} , for example, a constant in the range of 0.25 to 2.0 can be used.

(Rin)
Phosphorus elution is shown as in equation (83).

Here, PO ₄ (1) indicates the concentration of phosphoric acid in the aerobic layer. PO ₄ (2) indicates the concentration of phosphoric acid in the anaerobic layer.
As for phosphorus, the distribution coefficient π ₁ in the aerobic layer is expressed by the equation (84).

Here, _JP2 represents the decomposition phosphorus flux (unit: mg P / m ² / day) of the anaerobic layer.
Δπ _PO4,1 indicates a stepwise fractionation coefficient (unit: L / kg) of phosphoric acid in the aerobic layer. The value of Δπ _PO4,1 is set to 300, for example.
In addition, _πPO4,2 represents a step fraction coefficient (unit: L / kg) of phosphoric acid in the aerobic layer. The value of Δπ _PO4,1 is set to 100, for example.
[O2 (0)] _{crit, PO4} indicates the oxygen concentration in seawater (unit: mg / kg) that causes the fractionation coefficient in the aerobic layer to start decreasing stepwise. [O2 (0)] The values of _{crit and PO4} are set to 2.0, for example.
As for phosphorus, Expressions (51), (52), (53), and (54) are used as expressions indicating f _d1, f _p1, f _{d2, and} f _p2 , respectively.
(silica)
The elution of silica is shown as in formula (85).

Further, H ₂ (dP _Si / dt) is expressed as in Expression (86).

Here, P _Si indicates the concentration of suspended biological silica. S _Si indicates the generation rate of dissolved silica.
_JPSI is expressed as shown in Expression (87).

Also, _{k 2} is as shown in equation (88).

Further, π ₁ is expressed as in Expression (89).

Moreover, _kSi shows the dissolution rate constant (unit: d <^-1> ) of the particulate silica derived from algae. The value of k _Si is set to 0.50, for example.
Θ _Si3 indicates a temperature coefficient of silica oxidation. The value of θ _Si3 is set to 1.10, for example.
[Si] _Sat represents the saturation concentration of silica in the pore water (unit: mg Si / m ³ ). The value of [Si] _Sat is set to 40000, for example.

K _{M, PSi} represents a half-saturation constant (unit: mg Si / m ³ ) for the decomposition of particulate silica. The value of KM _{, PSi} is set to 5.0E + 07 (50000000), for example. Alternatively, mg Si / g may be used as a unit of _{KM and PSi} . In this case, the values of _{KM and PSi} are set to 100, for example.
Δπ _{Si, 1} indicates a stepwise fractionation coefficient (unit: L / kg) of silica in the aerobic layer. The value of Δπ _{Si, 1} is set to 10, for example.
Further, π _{Si, 2} represents a fractionation coefficient (unit: L / kg) of silica in the anaerobic layer. The value of π _{Si, 2} is set to 100, for example.
_JPSi represents the particulate silica flux from seawater to the bottom mud (unit: mg Si / m ² -d).
As for silica, Formula (53) is used as a formula indicating f _d2 .

<Water quality model>
In the water quality model, the change in the amount (concentration) of substances contained in water is calculated with a focus on the effects of biological activities on water quality. In the following, the carbon cycle related to diatoms will be described as an example.
FIG. 5 is an explanatory diagram showing an example of the circulation of substances related to diatoms simulated by the water quality model.
As shown in FIG. 5, the amount of diatom-related substances related to diatoms decreases as diatoms prey on bivalves and zooplankton. In addition, the amount of suspended organic matter increases due to the death of diatoms, and the amount of dissolved organic matter increases due to metabolism of diatoms. In addition, the amount of dissolved oxygen changes due to photosynthesis and respiration of diatoms. In addition, the amount of phosphate phosphorus, ammonia nitrogen, and silica is changed by respiration, killing, and photosynthesis of diatoms, and phosphate phosphorus is decomposed into nitric acid and nitrous acid.
The change in the amount of carbon related to diatoms is shown as in equation (90).

Here, Pdc indicates the carbon content of diatoms (unit: mgC / l). Pmd indicates the rate of increase in carbon content of diatoms (unit: day ⁻¹ ) under the optimum environment. The value of Pmd is set in advance based on, for example, observation or test. As the value of Pmd, for example, 2.25 can be used.
fd (T) represents the water temperature limiting function of diatoms. The water temperature limiting function fd (T) is a function that simulates the influence of the water temperature on the growth rate of diatoms. The growth of phytoplankton such as diatoms is activated until the environmental water temperature reaches the optimum water temperature, and the activity decreases when the temperature rises above the optimum water temperature. The relationship between this growth (activity) and water temperature change is similar to a normal distribution curve. Therefore, the water temperature limiting function fd (T) is set as shown in Expression (91).

Here, T represents the water temperature (unit: ° C.). Topt indicates the optimum water temperature (unit: ° C.). KTg1 represents a water temperature coefficient (unit: ° C. ⁻¹ ) related to growth below the optimum water temperature. KTg2 indicates a water temperature coefficient (unit: ° C ⁻¹ ) related to growth when the optimum water temperature is exceeded. The values of Topt, KTg1, and KTg2 are set in advance based on, for example, observation or test. For example, 20 (° C.) can be used as the value of Topt. For example, 0.0042 (/ ° C.) can be used as the value of KTg1. For example, 0.006 (/ ° C.) can be used as the value of KTg2.

  F (I) represents a light limiting function. The light limiting function f (I) is a function that simulates the influence of light on the growth rate of an organism. The light restriction related to the growth of phytoplankton is activated until the light amount reaches the optimum light amount, and the activity decreases when the light amount increases beyond the optimum light amount. Therefore, the light limiting function f (I) is set as shown in Equation (92).

  Here, ab is expressed as in Equation (93).

  Moreover, ax is shown like Formula (94).

Here, FD indicates the sunshine rate. The sunshine rate FD is set in advance based on, for example, weather data. For example, a constant in the range of 0.4 to 0.6 can be used as the value of FD.
Kess represents the total extinction coefficient of light (unit: m ⁻¹ ). The total extinction coefficient Kess of light is
As the existing amount of phytoplankton increases, the light absorption by the phytoplankton itself (light absorption by the cells) increases, and the light energy that can be used by the phytoplankton community in the deeper area is reduced. The total extinction coefficient Kess of light is expressed as shown in Equation (95).

Here, Keb represents an absorption coefficient in water (absorption coefficient of field) (unit: m ⁻¹ ). As the value of Keb, 0.7 can be used. Kechl represents an absorption coefficient (unit: m ² / mgchl) by phytoplankton (chlorophyll a amount). As the value of Kechl, 17 can be used. Pc represents the carbon content (unit: mgc / l) of macroalgae. Pdc indicates the carbon content of diatoms (unit: mgC / l) as described above. Pcc indicates the amount of carbon of dinoflagellates (unit: mgC / l). CChl indicates the ratio (unit: g C / mg chl) between the amount of carbon and the amount of chlorophyll a. The value of CCh1 is set in advance based on observations or tests, for example. As the value of CChl, for example, a value within the range of 5 to 100 can be used.
Δz represents the thickness of the calculation grid (that is, the distance between the grids in the vertical direction). ZD represents the distance (unit: m) from the water surface to the upper surface of the calculation grid.
Further, Is indicates an optimum light amount (unit: ly / day). Considering the time for the phytoplankton to adapt to the change in light, the optimum light amount Is is set as shown in Equation (96).

Here, Io indicates the amount of solar radiation on the water surface. I1 indicates the amount of solar radiation on the previous day (unit: ly / day). I2 indicates the amount of solar radiation (unit: ly / day) two days ago. The values of Io, I1, and I2 are set based on, for example, weather data.
Dopt indicates the maximum growth depth (unit: m) of phytoplankton. As the value of the maximum growth depth Dopt of phytoplankton, 1 can be used.
Fd (N) represents a nutrient limitation function. The nutrient limitation function fd (N) is a function that simulates the influence of the nutrient concentration on the growth rate of the organism, and is expressed as in Equation (97).

Here, KHn represents a half-saturation constant (unit: gN / m ³ ) related to nitrogen intake. As the value of KHn, for example, 0.01 can be used. KHp represents a half-saturation constant (unit: g P / m ³ ) related to phosphorus intake. As the value of KHp, for example, 0.001 can be used. NH ₄ and NO ₃ indicate the concentration of ammonia nitrogen and the concentration of nitrate nitrogen (both units are g N / m 3), respectively. PO ₄ represents the orthophosphate phosphorus concentration (unit: g P / m ³ ).

BMrd indicates the respiration rate (unit: day ⁻¹ ) at the optimum water temperature. As the value of the respiration rate Mrd at the optimum water temperature, for example, 0.01 can be used.
KTbm represents a water temperature coefficient (unit: ° C. ⁻¹ ) related to respiration. As a value of the water temperature coefficient KTbm related to respiration, 0.069 can be used. T represents the water temperature (unit: ° C.). The value of the water temperature T is set based on, for example, weather data. Topt indicates the optimum water temperature as described above.

BPrd indicates the killing rate (unit: ° C. ⁻¹ ) at the optimum water temperature. The value of the death rate BPrd at the optimal water temperature is set in advance based on, for example, observation or test. A constant in the range of 0.08 to 0.0215 can be used as the value of the killing rate BPrd at the optimum water temperature.
KTbp represents a water temperature coefficient (unit: ° C. ⁻¹ ) related to death. The value of the water temperature coefficient KTbp related to death is preset based on, for example, observation or test. For example, 0.069 can be used as the value of the water temperature coefficient KTbp related to death.

  Svpdc indicates the sedimentation rate of phytoplankton (unit: m / day). The value of the phytoplankton sedimentation rate Svpdc is set in advance based on, for example, observation or test. The value of the phytoplankton sedimentation rate Svpdc can be set to a constant value in the range of 0.02 to 0.5, for example.

ZPm represents the growth rate (unit: day ⁻¹ ) of zooplankton under the optimum environment. ZPm is represented by the formula (98) based on the existing amount Pcc of dinoflagellate, the existing amount Pdc of diatom, and the existing amount POC of particulate organic matter.

Here, ZPmax indicates the maximum value (unit: day ⁻¹ ) of the growth rate of zooplankton. The value of ZPmax can be obtained by experiment. Specifically, adult zooplankton females are collected and about 1000 eggs produced by the next morning in the laboratory are prepared as a group. While abundantly feeding the cultured phytoplankton as a bait, keep the temperature constant and measure the number of days required from hatching to adult growth to calculate the growth rate. Although the growth rate increases with increasing water temperature, it can be regarded as the maximum growth rate at this water temperature because the feed supply was obtained under abundant experimental conditions.
Prpfc indicates the phytoplankton intake rate. The value of the phytoplankton intake rate Prpfc is preset based on, for example, observations or tests. As the value of the phytoplankton intake rate Prpfc, for example, 0.5 can be used.
Prpfoc indicates the intake rate of suspended organic matter. The value of the uptake rate Prpfoc of the suspended organic matter is set in advance based on, for example, observation or test. For example, 0.5 can be used as the value of the suspended organic matter intake rate Prpfoc.
Zoomin indicates the minimum existing amount of zooplankton (unit: gC / m ³ ). When the zooplankton reaches zero, it cannot grow in the subsequent calculations, so the minimum existing amount is set. A value of 0.001 can be used as the value of the minimum zooplankton existing amount Zoomin.
Zsp2 indicates a half-saturation coefficient related to phytoplankton intake. The value of the half-saturation coefficient Zsp2 related to phytoplankton intake is set in advance based on, for example, observations or tests. As the value of the half-saturation coefficient Zsp2 related to phytoplankton intake, 0.3 can be used.
f (T) represents the water temperature limiting function of macroalgae. As described above, fd (T) represents a water temperature limiting function of diatoms. As f (T), fd (T) in formula (91) is used as f (T).
Zc represents the carbon content of zooplankton (unit: mgC / l).
Shc indicates the amount of intake by bivalves (unit: gC / m ² / day).

The product of “Pmd · fd (T) · f (I) · fd (N)” and “Pdc” on the right side of the formula (90) is the change in carbon amount due to diatom growth (here, the carbon amount related to diatoms). Change). The product of “BMrd · e ^{KTbm (T-Top)} ” and “Pdc” indicates a change in carbon amount due to respiration of diatoms. The product of “BPrd · e ^{KTbp (T-Top)} ” and “Pdc” indicates the change in carbon content due to the diatom death. The product of “Svpdc · (e / dz) and“ Pdc ”indicates the change in carbon content due to sedimentation of diatoms. “Zpm · f (T) · Zc” indicates a change in the amount of carbon caused by zooplankton preying on diatoms. “Shc” indicates the change in the carbon content due to the bivalve intake of diatoms as described above.

  The water quality model simulates the circulation of substances related to each of dinoflagellates, zooplankton and macroalgae, as in the case of diatoms. The water quality model simulates the circulation of organic carbon such as sedimentation of suspended organic carbon and decomposition of suspended organic carbon into dissolved organic carbon. The water quality model also simulates changes in the amount of organic matter related to land load. The water quality model simulates the circulation of phosphorus, the circulation of nitrogen, the circulation of dissolved oxygen, and the circulation of silica, as in the case of the circulation of carbon.

<Calculation of transparency>
The equation for calculating the transparency is shown as equation (99).

Here, the denominator “0.139 ד total SS concentration ”+ 0.019 ד chlorophyll a concentration ”+0.04” on the right side of the equation (99) represents the luminous flux extinction coefficient (K _d ). Therefore, the constant for each sea area is transparency × luminous flux extinction coefficient (Kd). By measuring the transparency and calculating the light flux extinction coefficient, a constant for each sea area can be obtained. For example, in the case of Tokyo Bay, the constant value for each sea area is set to 1.5. In the case of Ise Bay and the Seto Inland Sea, the constant value for each sea area is 1.6.
As shown in the equation (99), when calculating the transparency, it is necessary to obtain the total SS (Suspended Solid) concentration and the chlorophyll a concentration. The total SS concentration is calculated by a program calculation by dividing into two, an inorganic SS concentration and an organic SS concentration. The equation for obtaining the total SS concentration is shown as equation (100).

Here, as described above, POC, PON, and POP represent the particulate organic carbon concentration, the particulate organic nitrogen concentration, and the particulate organic phosphorus concentration, respectively.
As an inorganic SS, an inorganic SS from an inflow river is simulated. Specifically, the LQ equation for each river is calculated in advance at the time of low water and at the time of water discharge. The LQ equation is an equation showing the relationship between the river flow rate Q and the suspended sediment transport amount L.
The simulation calculation unit 191 calculates the floating sediment transport amount L by applying the river flow rate to the LQ equation, and calculates the amount of land-derived SS (inorganic SS flowing from the river) based on the floating sediment transport amount L. calculate. For example, 8% of the suspended sediment flowing in from the river is suspended sand and settles near the mouth of the river, so the simulation calculation unit 191 calculates that the remaining 92% flows into the simulation target water area.

  Further, the inorganic SS is divided into four types of less than 1 μm, 1 μm or more to less than 4 μm, 4 μm or more to less than 16 μm, or 16 μm or more according to the particle size. For example, the ratio of each section of less than 1 μm, 1 μm or more to less than 4 μm, or 4 μm or more to less than 16 μm is set to 2%, 8%, 50%, or 40%, respectively. The simulation calculation unit 191 calculates the inorganic SS amount for each section by multiplying the inorganic SS amount calculated based on the LQ equation by the ratio. Then, the simulation calculation unit 191 has a settling velocity (unit: m / day) for each section of less than 1 μm, 1 μm or more to less than 4 μm, 4 μm or more to less than 16 μm, or 16 μm or more (0), 0.28, Simulation is performed as 4.53 and 18.07.

The organic SS concentration can be obtained from POC, PON and POP. POC, PON, and POP are obtained for each layer and each mesh (for each lattice) in the water quality model. By summing up POC, PON and POP obtained for each layer and for each mesh, the organic SS concentration can be obtained.
Moreover, since all phytoplankton contains the chlorophyll a pigment | dye, the chlorophyll a density | concentration points out the quantity of the phytoplankton in seawater. For example, the amount of typical diatoms and dinoflagellates (carbon amount Pdc of diatoms and carbon amount Pcc of dinoflagellates) among phytoplankton can be calculated, and the total amount can be used as the chlorophyll a concentration.

Next, a program for executing the simulation will be described.
<Adoption of variable lattice>
FIG. 6 is an explanatory diagram illustrating an example of setting a variable grating. As shown in the figure, the grid of the portion to be obtained in detail is made fine, and a large grid is set in a place where importance is low. Thereby, a detailed calculation result can be obtained for a portion where the lattice is set finely. In addition, it is possible to reduce the amount of calculation for a place where a large lattice is set.

Here, the value of the distance (Δx and Δy) between the lattices is set as a variable of the two-dimensional array. In addition, variable lattices can be calculated by appropriately handling the inter-lattice distance used for calculation in the program.
FIG. 7 is an explanatory diagram illustrating an example of a program that uses the interstitial distance set as a variable of a two-dimensional array. This figure shows an example of a loop for obtaining the flux. The variable H1 indicates the interstitial distance (Δx) in the x direction. The variable H2 indicates the interstitial distance (Δy) in the y direction.
However, the locations using the interstitial distances set in the variables of the two-dimensional array are not limited to the locations shown in FIG.

<Improved stability calculation>
Rounding errors occur in computer calculations. If double precision floating point numbers are used, rounding errors can be reduced and the precision can be increased. However, if all items are calculated with double-precision floating-point numbers, computation time is required, such as taking 10 times or more time compared with the case of calculating with single-precision floating-point numbers, and a load is imposed on the computer.
Therefore, the water level, vertical eddy viscosity coefficient, and vertical eddy diffusion coefficient are expressed by double precision floating point numbers, and other data is expressed by single precision floating point numbers. As a result, the calculation can be performed with high accuracy, and the increase in the calculation time and the load on the computer can be suppressed as compared with the case where all the data is expressed in double precision floating point.

<Acceleration of calculation (1)>
In this model, most of the calculation time takes time in the loop processing of I, J, and K. The processing speed was improved by improving the loop. The improved loop, including the ecosystem, is over 100 and is treated in the same way. The original program is shown as an example below, and the subsequent improvement method is described.

FIG. 8 is an explanatory diagram showing an example of a program when the high-speed processing is not performed. The loops I, J, and K in the same figure indicate coordinates in the x direction, the y direction, and the z direction, respectively. In the program of FIG. 8, a loop is executed for the entire range of I and J.
FIG. 9 is an explanatory diagram illustrating an example of a program when the high-speed processing is performed. In the program of FIG. 9, the maximum values of I and J are stored in advance in INUM (n) and JNUM (n). Then, the loop is executed within the range of the stored maximum value. In this way, by suppressing the execution of the loop for the portion outside the range of the water area, for example, in the simulation of Tokyo Bay, the calculation time could be reduced to one third.

<Acceleration of calculation (2)>
The vertical loop was placed inside the horizontal loop to reduce the number of IF statement executions. Here, the vertical loop is a loop that sequentially processes lattice points in the vertical direction (z-axis direction). The horizontal loop is a loop that sequentially processes lattice points in the horizontal direction (x-axis direction and y-axis direction).
FIG. 10 is an explanatory diagram showing an example of a program in which the vertical loop is arranged inside the horizontal loop. In the program shown in the figure, the loop of variables I and J corresponds to an example of a horizontal loop, and the loop of variable K corresponds to an example of a vertical loop. The loop of variables I and J is executed by the loop of variable n starting at line L11 and ending at line L17. The values of variables I and J are set in rows L12 and L13, respectively. The loop for variable K begins at line L15 and ends at line L16.
The IF statement in the line L14 needs to be performed in a state where the values of the variables I and J are determined, and must be performed inside the loop of the variables I and J. On the other hand, the variable K does not appear in the IF statement in the line L12. For this reason, the IF statement in line L11 does not need to be inside the loop of variable K.

If the horizontal loop (variable I, J loop) is located inside the vertical loop (variable K loop), the IF statement will also be located inside the vertical loop. Must be repeated.
On the other hand, in the example of FIG. 10, the IF statement in the row L12 is located outside the vertical loop starting from the row L15. Thereby, in the example of FIG. 10, it is not necessary to repeatedly execute the IF statement every time the vertical loop is executed.

As described above, the IF sentence is located outside the vertical loop, so that the vertical loops can be parallelized relatively easily. Further, the CPU load can be reduced in that the number of execution times of the IF statement is reduced.
In the example of FIG. 10, the vertical loop is positioned inside the array declaration. Thus, optimization is performed so that vector calculation can be performed. Specifically, in the example of FIG. 9, AVD (MAXNUM, K) is used in the variable declaration at the top of the program. On the other hand, in the example of FIG. 10, the variable declaration at the top of the program is AVD (K, MAXNUM).
As shown in FIG. 10, the vertical loop is positioned inside the horizontal loop as shown in FIG. 10, and the vertical loop is also positioned inside the array declaration. In the experiment, the vertical loop is positioned outside the horizontal loop. The processing time was also reduced by two-fifths compared to the case where the vertical loop was located outside.

<Acceleration of calculation (3)>
Reduce processing time by parallelizing loops. In parallel, for example, OpenMP (Open Multi-Processing) technology is used.
FIG. 11 is an explanatory diagram showing an example of a program in which loops are parallelized. In the example shown in the figure, a loop of variables I and J is executed in a loop of variable n. Specifically, the loop of variable n starts at line L22 and ends at line L25. The values of variables I and J are set in rows L23 and L24 in this loop, respectively.
Further, parallelization is declared in line L21, and the loop of the variable n is executed in parallel. As a result, in the experiment, the speed could be increased by about 11 times using a 12-core CPU.

As described above, the model storage unit 181 is set for each grid set in the simulation target water area, and uses the projected area and volume of the pier included in the grid as parameters, and the governing equation for calculating the influence of the pier on the water flow. A model of water flow including The simulation processing unit 191 calculates a water flow in the simulation target water area based on the model stored in the model storage unit 181.
Thereby, in the water environment simulation apparatus 100, simulation can be performed with high accuracy even when a bridge is provided in the water area. In particular, in the water environment simulation apparatus 100, the projection area and volume can be used as parameters regardless of the shape of the pier, and therefore it is not necessary to customize the governing equation for each shape of the pier. In this respect, the burden on the user can be reduced.

The interstitial distance storage unit 182 stores the interstitial distance set for each position in the simulation target water area. The grid point coordinate calculation unit 192 calculates the coordinates of the grid point based on the inter-lattice distance stored in the inter-lattice distance storage unit 182.
Thereby, in the water area environment simulation apparatus 100, a lattice can be set flexibly according to the shape of the water area to be simulated.
As a result of the experiment, it was possible to stabilize the calculation by setting the change rate of the distance between adjacent lattices within ± 8 percent (%). Alternatively, the water environment simulation apparatus 100 or the user may set the lattice by calculating the inter-grid distance using a quadratic function with coordinates as variables (parameters). Thereby, a smooth water flow can be reproduced.

In addition, the data storage unit 183 stores data in which the water level, the vertical eddy viscosity coefficient, and the vertical eddy diffusion coefficient are expressed by numerical values having more digits than other numerical values. For example, the data storage unit 183 stores the water level, the vertical eddy viscosity coefficient, and the vertical eddy diffusion coefficient in double precision decimals, and stores other data in single precision decimals. And the simulation process part 191 calculates the water flow in the water area of a simulation object using the data which the data memory | storage part 183 memorize | stores.
Thus, by properly using the accuracy according to the type of data, the water area environment simulation device 100 can suppress an increase in necessary storage capacity and perform simulation with high accuracy.

The grid number storage unit 184 stores the number of grids for each row and column of the grid in the horizontal direction. Then, the loop processing unit 193 executes a loop in the row direction and a loop in the column direction within the range of the number of lattices stored in the lattice number storage unit 184. The simulation processing unit 191 calculates the water flow in the grid for each grid in the loop executed by the loop processing unit 193.
Thereby, in the water environment simulation apparatus 100, the number of executions of the loop can be reduced as compared with the case where the loop is executed a predetermined number of times regardless of the number of lattices. Thereby, in the water area environment simulation apparatus 100, the processing load of the water area environment simulation apparatus 100 can be reduced, and the simulation execution time can be shortened.

In addition, the loop processing unit performs grid row and column operations in the horizontal direction by parallel processing of the same loop.
Thereby, in the water environment simulation apparatus 100, the vertical loop can be arranged inside the horizontal loop, and the number of executions of the IF statement for determining the upper limit of the loop can be reduced. Thereby, in the water environment simulation apparatus 100, the parallelization rate of the loop can be increased, the processing load of the water area environment simulation apparatus 100 can be reduced, and the simulation execution time can be shortened.

In addition, the simulation processing unit 191 calculates water quality and pollution for each grid set in the simulation target water area, and calculates transparency for each grid point based on the obtained water quality and pollution.
Thus, according to the water area environment simulation device 100, the transparency in the water area to be simulated can be calculated, and the water quality of the water area can be evaluated based on this transparency.

DESCRIPTION OF SYMBOLS 100 Water environment simulation apparatus 110 Display part 120 Operation input part 130 Data input / output part 180 Storage part 181 Model storage part 182 Interstitial distance storage part 183 Data storage part 184 Grid number storage part 190 Control part 191 Simulation processing part 192 Grid point coordinate Value calculation unit 193 Loop processing unit

According to a first aspect of the present invention, water environment simulation apparatus, each grating set in simulated body of water, indicating the influence of the resistance by the pier the projected area and the volume of the pier which definitive within the lattice as a parameter A model storage unit that stores a model of a water flow including an equation including a term and an equation indicating an effect of turbulent energy on the mixing of the substance, and a simulation process for calculating a water flow in the simulation target water area based on the model A section.

The water area environment simulation device includes a data storage unit that stores data in which a water level, a vertical kinematic viscosity coefficient, and a vertical diffusion coefficient are represented by numerical values having more digits than other numerical values, and the simulation processing unit includes: You may make it calculate the water flow in the water area of the said simulation object using the data which the said data storage part memorize | stores.

The water area environment simulation device includes an interstitial distance storage unit that stores an interstitial distance set in a size corresponding to the position for each position in the simulation target water area, and a lattice based on the interstitial distance. A lattice point coordinate calculation unit that calculates the coordinates of the points, and the simulation processing unit calculates the lattice point coordinates by using the data stored in the data storage unit to simulate the water flow in the simulation target water area. It may be performed for each lattice point calculated by the unit.

The water environment simulation apparatus includes a grid number storage unit that stores the number of grids for each row and column of the grid in the horizontal direction , and the simulation processing unit corresponds to the rows and columns of the grid in the horizontal direction. Alternatively, the two-dimensional array may be replaced with a one-dimensional array based on the number of grids stored in the grid number storage unit , and the water flow in the grid may be calculated for each grid in a loop.

The water environment simulation apparatus may include a loop processing unit that executes the calculation of the rows and columns of the grid in the horizontal direction by parallel processing of the same loop.

The simulation processing unit may calculate transparency based on at least a total amount of carbon content of diatoms and carbon content of dinoflagellates .
The simulation processing unit may classify the suspended matters in the water into four stages according to the size of the particles, and calculate the concentration of each size of the particles.
The simulation processing unit may perform a simulation of decomposition of organic compounds, phosphorus compounds, and nitrogen compounds in the bottom mud and elution into water.

According to a second aspect of the present invention, water environment simulation method, each grating set in simulated body of water, indicating the influence of the resistance by the pier the projected area and the volume of the pier which definitive within the lattice as a parameter A water area environment simulation apparatus comprising a model storage unit that stores a model of a water flow including an expression including a term and an expression indicating an influence of turbulent energy on mixing of a substance, in the water area to be simulated based on the model A water flow state calculating step for calculating the water flow;

According to a third aspect of the present invention, the program for each grating set in simulated body of water, including a term indicating the influence of the resistance by the pier the projected area and the volume of the pier which definitive within the lattice as a parameter A water flow state for calculating a water flow in the simulation target water area based on the model in a computer having a model storage unit that stores a model of a water flow including a formula indicating an effect of turbulent energy on the mixing of substances A program for executing a calculation step.

Claims (8)

A model storage unit for storing a water flow model including a governing equation for calculating the influence of the pier on the water flow, using the projected area and volume of the pier included in the lattice as parameters, for each lattice set in the water area to be simulated; ,
A simulation processing unit for simulating water flow in the water area to be simulated based on the model;
A water environment simulation device. An interstitial distance storage unit that stores interstitial distances set for each position in the water area of the simulation target;
The water area environment simulation device according to claim 1, further comprising: a lattice point coordinate calculation unit that calculates coordinates of lattice points based on the inter-grid distance. A data storage unit for storing data in which the water level, the vertical eddy viscosity coefficient, and the vertical eddy diffusion coefficient are represented by numerical values having more digits than other numerical values,
The water area environment simulation apparatus according to claim 1, wherein the simulation processing unit calculates a water flow in the water area to be simulated using data stored in the data storage unit. A grid number storage unit that stores the number of grids for each row and column of the grid in the horizontal direction;
A loop processing unit for executing a loop in the row direction and a loop in the column direction within the range of the number of the lattices stored in the lattice number storage unit;
The water environment simulation device according to any one of claims 1 to 3, wherein the simulation processing unit calculates a water flow in the lattice for each lattice in a loop executed by the loop processing unit.   The water environment simulation device according to claim 4, wherein the loop processing unit executes the calculation of the rows and columns of the grid in the horizontal direction by parallel processing of the same loop.   The simulation processing unit calculates water quality and pollution for each grid set in the simulation target water area, and calculates transparency for each grid point based on the obtained water quality and pollution. The water environment simulation apparatus according to any one of 5.   For each grid set in the simulation target water area, a model storage unit for storing a water flow model including a governing equation for calculating the influence of the pier on the water flow using the projected area and volume of the pier included in the grid as parameters A water area environment simulation method, wherein a water area environment simulation apparatus includes a water flow state calculation step of calculating a water flow in the water area to be simulated based on the model. For each grid set in the water area to be simulated, a model storage unit for storing a water flow model including a governing equation for calculating the influence of the pier on the water flow using the projected area and volume of the pier included in the grid as parameters Computer
A program for executing a water flow state calculation step of calculating a water flow in the simulation target water area based on the model.