A novel model to simulate cloud dynamics with cellular automaton

Environmental Modelling and Software 122 (2019) 104537 Contents lists available at ScienceDirect Environmental Modelling and Software journal homepage: http://www.elsevier.com/locate/envsoft A novel model to simulate cloud dynamics with cellular automaton Alisson R. Silva, Adma R. Silva, Maury M. Gouv^ea Jr. * Polytechnic Institute, Departament of Electrical Engineering Pontifical Catholic, University of Minas Gerais, Belo Horizonte, Brazil A R T I C L E I N F O A B S T R A C T Keywords: Cellular automata Thermodynamic systems Cloud dynamics Dynamical systems Simulation In real system simulations, the application of cellular automata has been shown as an interesting option, because it can represent an emergent behavior and its implementation is simple. This paper presents a method for simulating thermodynamic systems, such as cloud dynamics, with cellular automata. In accordance with thermodynamic principles, this paper presents an isolated system model that describes temperature dynamics. The model uses the Von Neumann neighborhood of five cells, each with two possible states: the presence or absence of a cloud or a part of it. Our model uses three weather properties, as follows, condensed cloud water particles, temperature and outer winds. Four types of experiments were performed to validate the model proposed: (i) a warm body in the environment, (ii) a cloud dynamic simulation, (iii) a stability analysis, and (iv) a comparison between the dynamic of the proposed model and satellite images of clouds. 1. Introduction Simulations are studies of real systems using models that can describe and understand their behaviors and/or to evaluate strategies for their operations (Pedgen et al., 1990). Simulations can be performed in several areas, such as economics, biology, medicine, social sciences, and computation. There are some situations in which a simulation may be applied. For instance, when the plant of the process does not exist, when it is necessary to understand and analyze a real systems, when tests in a real system is not possible, when a real system can not be disturbed, or for training and education (Baik, 2005). Several advantages can be obtained with a simulation, in cost and time, because it can be repeated at lower cost or can represent a long time period in seconds. In order to simulate a real system, it is necessary to choose a simulation method (Fishman, 1996), (Gould and Tobochnik, 1988), (Rennard, 2002), (M. et al., 2010), (F. et al., 2010). Cellular automata (CAs) (Rennard, 2002) are discrete models on which many areas, such as computation, mathematics, physic, complexity science, and biology, are conducting research. CAs consist of a grid of cells, each with a finite number of states. For each cell, a set of cells, called neighborhood, is defined for the specified cell. At each iteration, a new state of each cell arises in accordance with the current state of the cell, the states of the cells in its neighborhood, and some fixed rules. Typically, the rule to update the state of a cell is the same for all cells and does not change over time, and is applied to the whole grid simultaneously. In real system simulations, the application of cellular automata has been shown as an interesting option, because it can represent an emergent behavior and its implementation is simple. CAs consist of a ndimensional grid of cells with the same behaviors described by a set of transition rules (Wolfram, 2002). CAs use a defined number of neighbor cells that interact with each other, creating a local interaction and then a global behavior. These interactions reflect the system dynamics based on the transitional rules. This paper presents a method for simulating thermodynamic systems, such as cloud dynamics, with cellular automata. A thermodynamic system is concerned with the flow and balance of energy. Three types of thermodynamic systems are distinguished depending on the types of interaction and energy exchange taking place between the system and its surrounding environment: (i) an isolated system is isolated in every way from its environment, and it does not exchange heat, work or matter with its environment; (ii) a closed system can exchange energy, heat and work, but not matter with its environment; and (iii) an open system exchanges energy and matter with its environment. A boundary that allows the exchange of matter is said to be permeable. In isolated systems, it is observed that as time passes, internal rearrangements decrease and stable conditions are reached. Properties, such as pressures and temperatures, tend to equalize, and matter arranges itself into one or a few homogeneous phases. A system in which all processes of change have ended is considered to be in a state of thermodynamic equilibrium. Clouds are formed from the condensation of water vapors present in * Corresponding author. E-mail address: [email protected] (M.M. Gouv^ ea). https://doi.org/10.1016/j.envsoft.2019.104537 Received 21 May 2018; Received in revised form 12 April 2019; Accepted 27 September 2019 Available online 2 October 2019 1364-8152/© 2019 Published by Elsevier Ltd. A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 the atmosphere. After it is formed, a cloud is moved by winds, and changes both its location and its properties, such as, temperature, pressure, density, and humidity. These properties strongly influence cloud dynamics. In accordance with thermodynamic principles, this paper presents an isolated system model that describes cloud dynamics. The model proposed uses the Von Neumann neighborhood of five cells, each with two possible states: the presence or absence of a cloud or a part of it. Our model uses three weather properties, as follows, condensed cloud water particles, temperature and outer winds. The transition rules are based on thermodynamic principles and weather concepts. Three types of experiments were performed: one with a warm body in the environment, other to analyze the system stability, and a cloud dynamics analysis. The latter experiments compare the evolution of the proposed model with satellite images of clouds. The rest of this paper is organized as follows. Section 2 presents important related works. Sections 3 and 4 introduce general concepts of cloud dynamics and cellular automata, respectively. Section 5 describes the proposed model of cloud dynamics model with cellular automata. Section 6 presents several experiments related to dynamics and stability analyses. Section 7 concludes the paper and suggests new implementations for future studies. Miyazaki et al. (2002) proposed a simulation technique, based a numerical solution of the partial differential equation of the atmospheric fluid model. Their model includes the phase transition and adiabatic cooling. The used the interaction of the vapor, the cloud, the temperature and the velocity field. In order to solve the partial differential equations they used the Boussinesq approximation, then assume that the air density is constant, so the atmospheric fluid is incompressible. To render the cloud, the model generate metaballs at the center of voxels. The images cloud are generate using the hardware-accelerated rendering method proposed by Dobashi et al. (2000). Their results showed that this model can simulate more realistic cumuliform clouds (cumulus and cumulonimbus) and clouds can be rendered that take the light scattering due to cloud particles into account. The model can simulate only cumuliform clouds. Harris (2003) developed a similar work to the work by Kajiya and Herzen and Overby et al. He developed a cloud dynamics simulation based on partial differential equations that model fluid flow, thermodynamics, and water condensation and evaporation, buoyant forces, and water phase transitions as well as various forces and other factors that influence these equations. But he also accounts for the negative buoyancy effects of condensed water mass and the positive buoyancy effects of water vapor. Other difference is that his model uses an exponential relationship between saturation and temperature. To solve the equations, he implemented the discrete form of these equations using programmable floating-point graphics hardware. All computation and rendering is performed on the GPU; the CPU provides only high-level control. The speed and parallelism of graphics hardware enables simulation of cloud dynamics in real time. His cloud model assumes that clouds exist alone, cannot represent clouds that are composed of tiny ice crystals or a mix of water and ice. The most important limitation of his cloud simulation system is the scale and detail that it can support. It is not currently possible to simulate a sky full of clouds. The domain size of most of the simulations that I have run is about 3–5 km at each dimension. These kind of approaches commonly needs a large amount of computation time because the equations produce a large of calculations and consequently a lot of data to process. In order to reduce the computational demands, other researchers have tried to simplify the cloud model using heuristic approach. Nagel and Raschke (1992) used a cellular automaton role-based model in order to simulate cloud dynamic. In the proposed model, the authors used a theoretical non mathematical description to develop a cellular automaton starting from the phenomenology of cloud formation. The authors used three logical variables: humidity, upward motion, and cloud water content. The state of these variables at the kþ1 instant is completely determined by the states of their neighbor cells at the k instant. Nagel and Rashke model uses only a small part of cloud physics process; nevertheless, the results showed important features of the fractal behavior of clouds and may be used as a starting point whenever a model of a cloud is needed, for instance, for problems of fractals analysis of satellite pictures (Rodriguez et al., 2012). Neyret (1997) applied some rule-base using a set of heuristics to simulate the dynamic effects of convective clouds. His method was fast to simulate, but the simulation effects were not realistic compared with physical method. Dobashi et al. (2000) extend four points to Nagel and Raschke model: extinction of clouds, wind effects, speeding up of the simulation and controlling cloud motion and they added a stochastic rule to evaporation and re-creation of clouds improving the model. Only cumulus clouds are simulated and rendered. These are more detailed and visually convincing. Miyazaki et al. (2002) used a simplified numerical model called Couple Map Lattice (CML) (Yanagita and Kaneko, 1997) to simulate cloud formation. Miyazaki et al. using rule based on atmospheric fluid dynamics considering an incompressible flow to solve Navier Stokes equation. To create rules, they have to be taken into account: viscosity 2. Related works Considering clouds and their simulation, there are various approaches that can be classified into two main groups: physically based and non-physically, based on heuristics, simulations. The physically based approach is closely connected to the simulation of fluid dynamics. A set of system of partial differential equations is required. The simulations based on heuristics use theoretical non mathematical description to create a role based systems. Some important works which used these two kind of approaches are presented as follow. Steiner (1973) developed a three-dimensional convection model. A set of partial differential atmospheric equations for motion, conservation, thermodynamic, and energy were solved using a staggered finite-difference technique. His model consider phase transition, but no-precipitation cloud. Because of the implicit nature of the set of equations for the buoyancy, pressure and saturation mixing ratio in deep moist convection, he used a shallow convection. A three-dimensional model was developed by Kajiya and Herzen (1984) to demonstrate a visual cloud simulation. Their focus was in rendering cloud using ray tracing algorithm. In order to generate cloud, they used a model that incorporate the equation of motion, continuity, condensation and evaporation. A set of partial differential equations, known as Navier-Stokes equations, define the model and these were solve numerically using incompressible flow, when it is assumed that the density of a fluid always remains constant over time, in this case fluid is the atmospheric air. However, this model does not include adiabatic cooling and the temperature lapse in simulation space, which is important for cumulus dynamics, also does not account for variation in pressure within the atmosphere and accounts for viscous effects, which are negligible for a standard atmospheric applications. Thus, the results is not so realistic. Overby et al. (2002) presented a similar work, but slightly more detailed. They combine the fluid solver with a model of the natural processes of cloud formation, including buoyancy, relative humidity, and condensation. To solve the Navier-Stokes equations, they used the stable fluid simulation algorithm of Stam (1999). The authors modified slightly processes to make the fluid model account for the variation in atmospheric pressure in the vertical dimension. They make two low-level modifications to their solver: the velocity field for the expansion of rising air and the conservation of momentum were modified. Although these effects tend to be minor (velocities in our typical simulation are modified by less than 5%), the effects are not negligible over the scale of the simulation. Furthermore, their model simulate more than one type of cloud. 2 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 and pressure effects, the advection of the state values by the fluid flow, diffusion of water vapor, thermal diffusion, thermal buoyancy, the phase transition from vapor to water. Their model could simulate various types of cloud such as cumulus, cumulonimbus, stratocumulus, altocumulus, and cirrocumulus. Also it obtained efficiently a realistic images of clouds by making use of graphics hardware. Guo et al. (2012) proposed a real-time dynamic cloud simulation approach based on multi-core and multi-thread. Differently from the works of Harris (2003), Kajiya and Herzen (1984), Overby et al. (2002), Dobashi et al. (2000), and Miyazaki et al. (2001) that focus on taking advantage of GPU computation and rendering ability, Guo et al. (2012) propose a cloud simulation method that both take advantage multi-core CPU and modern GPU’s parallel ability to achieve better computational performance and realistic dynamic cloud simulation. They used the Dobashi et al. (2000) idea in dynamic cloud modeling, but they introduced horizontal single-direction wind function to simulate cloud movement, to simply computation and enhance real-time performance. In order to render the clouds, they used Harris method. They adopted multi-thread to implement cloud modeling and illumination modeling to run on CPU while the illuminate color of clouds is rendered on GPU. The results showed a better performance than other similar work, and the better scalable ability on different multi-core platform. Qiu et al. (2013) proposed an effective method for simulating 3D cloud. They used CML model proposed by Miyazaki et al. (2001) to simulate cloud generate. They presented a new approach to simulate light scattering in clouds based on spherical harmonics and spherical harmonic coefficients that represent incident-light distribution in order to improve rendering a frequency domain volume-rendering algorithm combine with spherical harmonics was implemented. Their method was divided into two phases: off-line and real-time. To generate cloud data, pre-compute incident light and transform volume data to the representation of frequency domain is solved off-line. The work of real-time rendering processing is to carry out frequency domain volume rendering. The calculation process is divided in two phases: CPU processing and GPU processing. On CPU is generate a Cloud Data that is translated into 2D texture slices. In the follow calculation process, the spherical harmonic coefficients are updated according to iterative calculation results. On GPU, the rendering process is executed. The spherical harmonic coefficients are calculated according to information received from CPU. The results showed that their method facilitates computing efficiency while maintained a realistic visual quality. This method not considering flying through or inside the clouds, then is not very suitable for animation of clouds. dynamic, thermodynamic, and cloud physical. The equations of the cloud dynamics are governed by Newton laws of motions and continuity equations applied to the air, an equation for temperature, and conservation equations for water. The basic elements necessary to simulate clouds are air velocity, u ¼ ðux ; uy ; uz Þ, atmosphere pressure, p, temperature, T, water vapor, w, and condensed cloud water, μ. These water vapor and condensed cloud water variables are mixed ratios, i.e., the mass of vapor or liquid water per unit mass of air. In order to simulate the cloud dynamics in two dimension of a layer of the atmosphere as a function of the environment temperature, it is required four partial differential equations (Rogers and Yau, 1988): one of them for horizontal air speed, 3. Cloud dynamics 4. Cellular automata Clouds are formed from the condensation of water existing in the humid air in the atmosphere. The elevation of air is the key process in the production of clouds, because when it rises and comes into contact with low temperatures, cold air makes it possible for clouds to form. This elevation can be produced by convection, convergence of air streams, topographical elevation or frontal lifting (Vianello, 2000). Clouds may be in a liquid or solid state, or may be a mixed composition of water and ice. The composition of a cloud depends on its altitude. After having formed, clouds are moved by winds in all directions. When a cloud is moved in a vertical direction, its altitude changes as do its properties, such as temperature, pressure, kinetic energy, density, and humidity. On rising, there is a cooling of condensed cloud water particles that may become partially or completely frozen. On the other hand, when a cloud goes to a lower altitude, it goes to a higher temperature environment; therefore, precipitation may arise and spread the cloud. The dynamics, growth, motion and dissipation of clouds are complex. Thus, it is important to understand these dynamics in order to allow an efficient implementation of the real system (Andrews, 2000). Three type of processes can be accounted for a cloud model, as follows, A cellular automaton is formally defined as a discrete mathematical model, implemented in computers, automated by deterministic rules, and its conduct of an element within a homogeneous set will be based both on the state of its own attributes and those of the neighboring elements (Wolfram, 2002). A CA is characterized by its cell space and its transition rule. The cell space is a lattice of N identical cells arranged in a d-dimensional grid, each with an identical pattern of local connections to other cells. When we consider the lattice is of finite length, boundary conditions are applied resulting in a circular lattice. A transition rule provides the next state for each cell, as a function of the configuration of its current neighborhood. At each step of time, every cell of the lattice updates its states according to this rule. As to the CA-dimensional rule contained in each cell, it is essentially a finite state machine, usually specified in the form of a table of rules. These are called elementary cellular automata. The neighbors of a cell are adjacent cells, or cells on the right and left. Thus a cell is connected to air local neighbors (cells) where r is related to the radius, so that each cell has a neighborhood of 2r þ 1. A neighborhood is made up of three cells, so there are 23 ¼ 8 possible patterns for a neighborhood. There are ∂ux ¼ ∂t ux ∂ux ∂x uy ∂ux ∂y 1 ∂p þ Fu x ; ρ0 ∂x (1) one for vertical air speed, ∂uy ¼ ∂t ux ∂uy ∂x uy ∂uy ∂y 1 ∂p þ Fuy ; ρ0 ∂ y (2) one for continuity equation, 0¼ � ∂ ∂ ðρ u Þ þ ρ u ; ∂x 0 x ∂y 0 y (3) one for temperature, ∂T ¼ ∂t ux ∂T ∂x uy ∂T þ FT þ φT ; ∂y (4) one for water vapor, ∂w ¼ ∂t ux ∂w ∂x uy ∂w þ Fw þ φ w ; ∂y (5) and one for cloud water, ∂μ ¼ ∂t ux ∂μ ∂x uy ∂μ þ Fμ þ φμ : ∂y (6) The Fi and φj variables are buoyant force and entropy on i ¼ fux ; uy ; T; w; μg and j ¼ fT; w; μg variables, respectively. In the continuity Equation (3), the local rate of change of density has set to zero. This is called the anelastic assumption, which assumes that ρ0 is everywhere constant except when buoyancy is calculated – hence the continuity equation is eliminated. 3 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 therefore 28 ¼ 256 possible rules. Wolfram (2002) proposed a numbering scheme for the elementary CAs, in which the output bits are ordered alphabetically, as in the transition rule, and are read from right to left to form a base number in decimal notation between 0 and 255. For CAs, dimensional cells are arranged in a two-dimensional space (represented in the form of a grid), the neighborhood most widely used are the Von Neumann neighborhood, consisting of 5 cells, i.e., a central cell and 4 neighbors: up, down, left, and right) and the Moore neighborhood, consisting of 9 cells, i.e., the Von Neumann neighborhood of more cells in the diagonal. Cellular automata are used in simulation and emulation of real systems (Rennard, 2002), such as: � Simulation of bacterial or viral behavior, crystal growth, coral, rocks and other natural elements, behavior of gases, spread of fires, population development, economic, behavior of land, rivers and topographies, and forecast of plant growth; � Video: generating random pictures, image filters and distortions; � Music: melody-generating digital noise and sound; � Mathematics: alternative to replacement differential equations by difference equations; � Computer: random number generation, cryptography, and conceptual design of parallel computations mass; and � 3D animation: particle simulation and generating textures. Fig. 2. Vectorial field representing a wind flow. 5.1. Dew point temperature The dew point is defined as the temperature in which the air must be cooled to the condensation of water, i.e., the air becomes saturated with water vapor. At the temperature of the dew point, the amount of water vapor in the air is currently at its maximum concentration. The dew point temperature, Td , can be calculated as proposed by Vianello (2000) 5. Cloud dynamics with cellular automaton In our model, the interactions between cloud dynamics and microphysics is through two-dimension way. The cloud water is not present as described by Equation (5), but from the condensation of water vapors present in the atmosphere. Thus, a cloud appears when the temperature decreases and water vapor condenses to form cloud water. The limit used to condense water vapors is the dew point temperature, Td . This is the threshold in which the amount of water vapor in the atmosphere is currently at its maximum concentration (Vianello, 2000). The proposed model simulates clouds into a cellular automaton grid. The cloud has an initial size that may be modified by weather events, such as, temperature and wind. The latter is represented, into the grid, by a vectorial field. A cell, which represents a micro-region in the atmosphere, has three attributes: the temperature, the vertical and horizontal components of the wind vector. The following sections explain these parameters and the cloud dynamics. The model uses the von Neuman neighborhood with five cells, as shown in Fig. 1. Every cell has a ðx; yÞ coordinate, i.e., its position in the simulated area. Thus, xl ;xr ;ya , and yb are coordinates of the i-th cell neighbors, and Δx and Δy are the lengths between two horizontal and vertical neighbors, respectively. The Δx and Δy values depend on the number of cells used into the grid for an area. The higher is the number of cells into a grid, the lower Δx and Δy are. The next subsections describe the dynamic of the proposed model. Td ¼ 186:4905 237:3 loge ; loge 8:2859 (7) where e ¼ es A p ðT Tu Þ (8) is the real vapor pressure, es is the water vapor saturation pressure, A is the psychometrics constant, p is local atmospheric pressure, and (T Tu ) is the psychometrics depression. In our model, a cloud appears into the grid of the cellular automaton when the atmospheric temperature belongs to ½Td δ; Td þδ� interval, where δ is a constant. 5.2. Wind The wind flow is described by a vectorial field, where each vector is implemented into a cell and described by horizontal, ux , and vertical, uy , components. Fig. 2 shows a 5 � 5 cellular automata grid covered by a wind flow with a general direction from the south-west to the northeast. A general wind flow direction is established, but each vector component has a small disturbance to guarantee randomness, as can be seen in some vectors which have a bit different direction from the general one. The dynamics of the wind flow is defined by Equations (1) and (2). In discrete form, the partial derivatives of the wind speed and pressure with respect to the variables x and y, for the i-th cell at a given iteration, are approximated, respectively, by ∂ux Δuxi ; � ∂x Δxi (9) ∂uy Δuyi ; � ∂y Δyi (10) ∂p Δpi ; � ∂x Δxi (11) and ∂p Δpi ; � ∂y Δyi Fig. 1. The von Neumann Neighborhood. 4 (12) A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 pseudo-code of our model is presented as follows: where Δuxi ¼ ul ur , Δpi ¼ pl pr , and Δxi ¼ xl xr for Equation (9), and Δuyi ¼ ua ub , Δpi ¼ pa pb , and Δyi ¼ ya yb for Equation (10). Thus, the linearization of Equations (1) and (2) by Euler’s method provides the following difference equations � � Δuxi Δuxi 1 Δpi uxi ðt þ ΔtÞ ¼ uxi ðkÞ þ uxi þ Fuxi Δt; uyi (13) Δxi Δyi ρ0 Δxi for the horizontal wind speed, and � Δuyi Δuyi uyi ðt þ ΔtÞ ¼ uyi ðkÞ þ uyi uxi Δxi Δyi 1 Δpi þ Fuyi ρ0 Δyi � Δt; For a given altitude, initialize ρ0 , Td , and Ti and pi , 8i; Initialize Fi , i ¼ fux ; uy ; Tg, and φT ; Introduce a vectorial field of wind into the grid; Introduce a cloud into the grid; Update cell temperatures by interaction to its neighborhoods and the wind flow, according to Equation (19); 6. Execute the wind flow dynamics according to Equations (13) and (14); 7. Go back to Step 5 while a stop criterion is not satisfied. 1. 2. 3. 4. 5. (14) A microregion of the atmosphere will have a cloud if the i-th cell temperature, Ti , belongs to ½Td δ; Td þδ� dew point interval. for the vertical wind speed. 5.3. Temperature 6. Simulations and result analyses The i-th cell temperature, Ti ðtÞ, defines the value of this variable into this micro-region of the space. The value of Ti ðtÞ defines the appearance of cloud into the cell. Each cell has a cloud, or part of it, if Ti ðtÞ belongs to ½Td δ; Td þδ� interval; otherwise, there is no cloud into the cell. The cell temperature behavior is derived from the Newton’s law of cooling, in which heat transfer between the i-th cell and its neighborhood is provided by dTi ¼ dt q ðTi Tn Þ; This section presents four different experiments to validate the proposed model. The first experiment shows the performance of the model describing an isolated system, the second one shows two experiments with a 100 � 100 grid in order to present the cloud dynamics, the third experiment presents a stability analysis of the proposed model, and the last experiment compares the evolution of the proposed model with satellite images of clouds. (15) 6.1. Thermodynamic model with cellular automata where Tn is the neighborhood cell temperature and q is a positive constant. Since the i-th cell has four neighborhood, the model approximates dTi =dt as an average of Equation (15) for all neighborhood, as follows � � dTi 1 � q ½ðTi ðkÞ Tr Þ þ ðTi ðkÞ Tl Þ þ ðTi ðtÞ Ta Þ þ ðTi ðtÞ Tb Þ� No 1 dt i X q h � Tj ðtÞ ; 4 Ti ðtÞ No 1 j This section presents an experiment to describe an isolated system with cellular automata. In isolated systems, as the time passes, internal rearrangements decrease and stable conditions are reached. Properties, e.g., pressures and temperatures, tend to equalize, such that the processes of change come to an end and the system reaches the state of thermodynamic equilibrium. Two simulations of isolated thermodynamic systems are presented in 30 � 30 and 50 � 50 grids. The model was implemented with the Visual Studio 2010, using Cþþ language. For these simulations, the temperature interval was [0, 50] degrees, the initial temperature was 0 � C, and Δt ¼ 0:01. Fig. 3(c) shows the temperature intervals and their respective colors. Fig. 3 shows six iterations of each simulation. In the simulation of the 30 � 30 grid, Fig. 3(a), a warm body of 50 � C and 50%-grid area was initialized into the center of the grid. The heat of the warm body spread quickly – see iteration 4. In iteration 27, the extreme area of the grid had been warmed, by heat from the warm body, to a temperature of between 1 and 10 � C. Iteration 48 shows that the center of the warm body started to cool, because its heat spread throughout the grid. Iterations 112 and 369 show that the heat continued to spread, thus providing the temperature equalization. Iteration 472 shows the moment at which the grid temperature was totally equalized between 1 and 10 � C. In the simulation of the 50 � 50 grid, Fig. 3(b), a warm body of 50 � C and 80%-grid area was initialized into the center of the grid. The system behavior was similar to that of the first simulation, 30 � 30 grid. The heat spread quickly throughout the whole grid. After iteration 46, the heat of the center of the grid started to decrease. Iterations 112 and 134 show the beginning of the temperature equalization throughout the whole grid. Iteration 485 shows the grid with its temperature totally equalized between 20 and 30 � C, higher than that of the first simulation because of the larger area of the warm body (80%-grid area). These simulations showed that a cellular automaton model can simulate a thermodynamic system. In both of them, the heat in the center of grid spread throughout the whole grid until the thermodynamic equilibrium. (16) where j ¼ fl;r;a;bg, No is the number of cells per neighborhood, Ti ðkÞ is the temperature into i-th cell at iteration k, and the subscripts l, r, a, b mean the neighbors on the left, on the right, above, and below neighbors of an i-th cell. The cell temperature also interacts to the wind flow, and its dynamics is approximated in discrete time for each cell according to Equation (4). The partial derivatives of T with respect to x and y, for the i-th cell, are approximated by ∂T ΔTi � ∂x Δxi (17) and ∂T ΔTi ; � ∂y Δyi (18) respectively, where ΔTi ¼ Tl Tr and Δxi ¼ xl xr for Equation (17), and ΔTi ¼ Ta Tb and Δyi ¼ ya yb for Equation (18). Thus, the i-th cell temperature is updated as follows � � ΔTi ΔTi dTi Ti ðt þ ΔtÞ ¼ Ti ðkÞ þ Δt: (19) uyi þ FT þ φT þ uxi Δxi Δyi dt where dTi =dt comes from Equation (16). 5.4. The model dynamics The model dynamics is based on a discrete and iterative system. All grid cell transitions are based on Equations (13), (14) and (19) that change each cell state locally, providing the global effect in the grid. The 6.2. Cloud dynamics This subsection presents two simulations with a 100 � 100 grid, the 5 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 Fig. 3. Isolated thermodynamic system simulations. Fig. 4. From initial state to final one of the centered cloud dynamics simulation. 6 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 Fig. 5. From initial state to final one of the top right cloud dynamics simulation. Fig. 6. Atmospheric temperature evolution in cloud dynamics simulations. first one with the cloud in the center of the grid and other with the cloud in the top right of the grid. The grid altitude was 5000 m, the current atmospheric temperature was 3o C, the dew point Td ¼ 2 oC, Δt ¼ 0:01, and a fixed atmospheric pressure of 700 hPa. The cloud experiment ran 5000 iterations. Figs. 4(a) and 5(a) show the initial state of the two simulations. The cloud dynamics of our model showed an expected behavior regarding some thermodynamics concepts, because in all experiments the whole grid reached thermal equilibrium resulting in total cloud dissipation to which the actions of wind also contributed, as can be seen in Figs. 4 and 5. These figures show that the clouds were dispersed and moved due to the thermal equilibrium and wind interactions effects. In Figs. 4(d) and 5(d), the two clouds are smaller than their initial positions, Figs. 4(a) and 5(a), by the dissipation, and far from it, due to the wind interactions. In a grid where a given temperature prevails over almost of its whole area, the environmental temperature set to 3 � C, is expected at infinite time into an undisturbed environment such that it reaches thermal equilibrium close to the initial environmental temperature. Fig. 6 shows both the mean and standard deviation (SD) of the atmospheric temperature in the two experiments. Fig. 6(a) shows the mean temperatures of two experiments, which started from about 2.4 � C and finished to about 3 � C, almost the initial atmospheric temperature, i.e., both of them tend to thermal equilibrium. Fig. 6(b) presents SD of the two temperatures, which passed from 1.5 to about zero, i.e., in the two experiments the grid tended to a same temperature close to 3 � C. These effects resulting from physical phenomena were more prevalent than those of the simulations with smaller grids. These results may 7 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 be expected because of these smaller areas, thermodynamic equilibrium tends to be reached faster than in those of the larger grids, in addition to which the probability of a wind reaching the cloud is greater. unstable. On the other side, the precision of the system increases as Δt decreases; however, the period of simulation also decreases, becoming the simulation not feasible for real time. In this stability analysis, we assessed the largest Δt in which the system remain stable. Thus, the larger a stable Δt, the longer the period of simulation. In order to use an as precise as possible Δt, a deviation analysis between outputs from different simulations was performed. Fig. 7 presents the mean temperatures in the 100 � 100 matrix for different Δt. In Fig. 7(a), it is observed an unstable system, since Tm increases linearly to 1.8 � 1024 at iteration 10. Several iterations after that, Tm overflow (∞). Thus, it is expected that Tm →∞ as t→∞. In Fig. 7 (b), Tm keeps constant around the initial atmospheric temperature until the first 40 iterations; after that, Tm decreases linearly to ∞. Thus, it is expected that Tm →∞ as t→∞. In Fig. 7(c) and (d), we see a stable system that reaches the thermal equilibrium, with final Tm close to the initial atmospheric temperature. In Fig. 7(c), the mean temperature reaches the atmospheric temperature at 600 iterations. In Fig. 7(d), Tm stabilized at 5000 iterations. Fig. 8 presents the mean temperatures of the 2000 � 2000 matrix for different Δt. In Fig. 8(a) and (b), the system had similar unstable behavior of those presented in Fig. 7(a) and (b). Fig. 8(a) shows that Tm performed an inverse behavior with respect to that presented in Fig. 7 (a), i.e. in the former, Tm decreases linearly to ∞ (overflow) instead of to increase to ∞, in the latter. Similarly, in Fig. 8(b), Tm has an inverse behavior with respect to that in Fig. 7(b). In Fig. 8(b), the temperature kept constant until 40 iterations; after that, Tm increased linearly to ∞ (overflow). In Fig. 8(c), Tm stabilized with 8200 iterations and kept constant until the end of the experiment. Thus, it is possible to assume that the system reaches the thermal equilibrium with Δt ¼ 1 ms. In 6.3. Stability analysis This section presents the computational experiments on stability analysis of the proposed model. The stability is defined with respect to the fixed point of the system. When the system reaches a fixed point, T � , its state does not change as the time goes to infinity, i.e., TðtÞ→ T � as t→ ∞. In this work experiments, the stability was analyzed empirically, that is, we assumed that the system is stable whenever the mean temperature of the matrix does not change until the end of the simulation (20,000 iterations) and its value, TðtÞ, is close to the initial atmospheric temperature. Thus, we assume that the system reached the thermal equilibrium. Some parameters of the system were fixed in all experiments, and their values were defined as follow: q ¼ 0,0243 w/m 1 k 1, cell area was 0.01 m2, altitude was 5 km, atmosphere temperature was 3 � C, Td ¼ -8.4 � C, atmospheric pressure was 700 hPa, and wind speed was 60 km/ h. In order to analyze the system stability, several matrix dimensions and time step size, Δt, were used, totaling 16 experiments. The following matrix dimensions 100 � 100, 2000 � 2000, and 5000 � 5000, and Δt ¼ { 1 s, 0.1 s, 1 ms, 0.1 ms } were used. In discrete models, Δt determines the system precision and stability, and the real period of simulation. A large Δt provides a long period of simulation. Nevertheless, large Δt increases the deviation between the real and simulated system output. Moreover, the system may become Fig. 7. Mean temperature of the experiments with 100 � 100 matrix: (a) Δt ¼ 1 s (b) Δt ¼ 0.1 s (c) Δt ¼ 1 ms (d) Δt ¼ 0.1 ms 8 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 Fig. 8. Mean temperature of the experiments with 2000 � 2000 matrix: (a) Δt ¼ 1 s (b) Δt ¼ 0.1 s (c) Δt ¼ 1 ms (d) Δt ¼ 0.1 ms Fig. 8(d), it is not possible to reach a conclusion. Nevertheless, likely the system will stabilize, since in 20,000 iterations Tm decreased a little due to the disturbance caused by the cloud (colder than the initial atmospheric temperature). By the reason of the matrix size, 20,000 iterations were not enough to stabilize the system. Fig. 9 shows the mean temperatures of the 5000 � 5000 matrix for several Δt. Fig. 9(a) shows the experiment where the behavior of the system was unstable with Tm higher than 3 � 1021 in only 10 iterations. Such as in the last experiments, we assume that Tm →∞ as t→ ∞. In the experiment shown in Fig. 9(b), the system had similar behavior to that in 100 � 100 matrix with Δt ¼ 0:1 s, Fig. 7(b). In Fig. 9(c), Tm stabilized about 19,000 iterations. Thus, the system stabilized at a slow rate likely due to the matrix dimension. In the experiment shown in Fig. 9(d), the system had a similar behavior to that presented in Fig. 8(d). Such as in that experiment, we assume that the matrix dimension and low Δt did not permitted that 20,000 iterations were enough to stabilize the system. In Fig. 7, Δt equals 1 ms and 0.1 ms provided a stable behavior. In Figs. 8 and 9, only Δt ¼ 1 ms provided a stable behavior. In all experiments, Δt equals 1 s and 0.1 s provided an unstable behavior; then, these two value can be considered not suitable for the proposed model. The system behaved stable in all experiments with Δt ¼ 1 ms until 20,000 iterations. A deviation analysis of the mean temperatures with Δt equals to 1 ms and 0.1 ms was performed with 100 � 100 matrix in order to check the evolution of the proposed model. As shown in Fig. 10, until iteration 100 the values of the mean temperatures were similar; after that, there was a deviation equals to 0.002 until iteration 280; then, the deviation increased to 0.0063 until iteration 580. Finally, the deviation decreased to zero until the end of the simulation. This error can be considered acceptable with respect to the magnitude of the temperature. Table 1 shows the means and standard deviations of the temperatures with Δt ¼ 1 s for all matrix dimensions. The initial mean temperatures decreased as the matrix dimensions increased due to the cloud area with respect to the matrix dimension, i.e. the larger is the matrix dimension, the smaller is the disturb caused by the cloud in the environment. The initial standard deviations also show this relation. For instance, in 100 � 100 matrix, the rate cloud area per matrix area is larger than that in 5000 � 5000 matrix; thus, in the former, the cloud provides more disturbance. It can be seen that the final mean temperatures and standard deviations provided nan (not a number) value, i.e. these variables went to �∞. These results show that the system did not reach the thermal equilibrium, in accordance with the analysis on Figs. 7 to 9. Table 2 shows the means and standard deviations of the temperatures with Δt ¼ 0.1 s, in which they had similar behavior to those with Δt ¼ 1 s, Table 1. Table 2 also shows an unstable system, where the final Tm and standard deviations provided nan, i.e., they tended to ∞. Table 3 shows the means and standard deviations of the temperatures with Δt ¼ 1 ms. In all experiments, Tm and standard deviations tended to the initial atmospheric temperature and zero, respectively. These data support our analysis based on Figs. 7(c), 8(c) and 9(c). 9 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 Fig. 9. Mean temperature of the experiments with 5000 � 5000 matrix: (a) Δt ¼ 1 s (b) Δt ¼ 0.1 s (c) Δt ¼ 1 ms (d) Δt ¼ 0.1 ms Table 1 Mean temperatures and standard deviations with Δt ¼ 1 s. Dimension 100 � 100 2000 � 2000 5000 � 5000 Mean Temperature Standard Deviation Initial Final Initial Final Time (s) –nan –nan –nan 1.7754 0.5671 0.3225 –nan –nan –nan 16.857 5486.146 35,540.965 3.6718 3.0605 3.0194 Table 2 Mean temperatures and standard deviations with Δt ¼ 0.1 s Dimension 100 � 100 2000 � 2000 5000 � 5000 Mean Temperature Standard Deviation Initial Final Initial Final Time (s) –nan –nan –nan 1.7754 0.5671 0.3225 –nan –nan –nan 16.721 5487.250 35,549.855 3.6718 3.0605 3.0194 Table 3 Mean temperatures and standard deviations with Δt ¼ 1 ms Fig. 10. Deviation of the mean temperatures in 100 � 100 matrix provided by Δt ¼ {1, 0.1} ms. Dimension Finally, Table 4 shows the means and standard deviations of the temperatures with Δt ¼ 0.1 ms. The final mean temperature and standard deviation in 100 � 100 matrix tended to the initial atmospheric temperature and zero, respectively. In this experiment, we can assume that the system reached the thermal equilibrium. In experiments with 2000 � 2000 and 5000 � 5000 matrices, the mean temperatures did not 100 � 100 2000 � 2000 5000 � 5000 10 Mean Temperature Standard Deviation Initial Initial Final Time (s) 1.7754 0.5671 0.3225 0.00050 0.00034 0.00028 17.251 5516.543 35,653.547 3.6718 3.0605 3.0194 Final 3.00004 3.00003 3.00004 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 Fig. 14(b). This occurred likely because of the convection between the atmosphere layer where the simulation takes place and the above or below one. Cloud (C), in Fig. 14(a), begins to dissolve due to heat exchange. Cloud (B), in Fig. 14(a), is seen to have a similar shape to its equivalent one in Fig. 14(b), except that it is less displaced by the wind. Finally, cloud (R) spreads in the atmosphere more than its equivalent in Fig. 14(b); however, their positions and shapes are similar. Table 4 Mean temperatures and standard deviations with Δt ¼ 0.1 ms Dimension 100 � 100 2000 � 2000 5000 � 5000 Mean Temperature Standard Deviation Initial Initial Final Time (s) 1.7754 0.5671 0.3225 0.0005 0.5207 0.3041 16.977 5497.412 35,583.871 3.6718 3.0605 3.0194 Final 3.0001 3.0606 3.0195 7. Conclusion change significantly along the simulations. In these experiments, the standard deviations showed a small decrease along the time. There are two reasons for that: the very small Δt, which evolves a small period of time during the simulation; and the matrix dimension, which needs many iterations to propagate the temperature disturbance. The number of time step size, Δt, used in the experiments provided the following conclusions: (i) values of Δt equal or larger than 0.1 s likely provide an unstable system, such as those used in our experiments – Δt ¼ { 0.1 s, 1 s }; and (ii) values of Δt smaller than 0.1 ms will provide a small evolution of the system for a long time of processing. Simulations with small Δt will not produce significant variations on the state of the system, especially in large matrix, such as those used in the experiments with Δt ¼ 0.1 ms for 2000 � 2000 and 5000 � 5000 matrices. The simulations showed that Δt ¼ 1 ms provides an useful step of time for the proposed model. This Δt supplies a stable system, with lead the area of simulation to the thermal equilibrium. This article set out to construct a model to simulate thermodynamic systems using cellular automata. Two types of models were presented, an isolated thermodynamic system and a cloud dynamics model. The former showed that a cellular automaton can simulate a thermodynamic system. This first model was the basis for the second one, the cloud dynamics model. The second model used a limited representation considering the variables that represent the dynamics of a real cloud. The atmospheric temperature and wind flow in the two-dimensional grid were included. A two-dimensional CA with a Von Neumann neighborhood of 5 cells was used. The transition rules were defined based on the thermodynamic principle, that defines the thermal equilibrium, and the wind interaction. The validation of the second model was made a simulation with a 100 � 100 cellular automaton grid. In this preliminary study, we did not compare our model with other ones for cloud dynamics. Thus, we only conducted visual validation and graphical analyses, considering basic thermodynamic principles. The simulations showed that our proposed model presented a satisfactory behavior, considering some thermodynamic principles. The wind actions also were considered coherent, because they moved the clouds until their full dissipation. In the experiment, the clouds, with heterogeneous temperature randomly initialized, tended to converge to a uniform temperature, reaching the thermal equilibrium. Another observed behavior was the thermal equilibrium between the cloud and grid, which always resulted in cloud dissipation. The wind actions also contributed to convergence of cloud temperature to environmental temperature, because they spread the clouds, accelerating the thermal equilibrium. Moreover, the proposed model was compared with realist data, using satellite images. The results showed that the dynamic of the proposed model was similar to that of the real satellite images. As a follow-up to this study, other variables will be added, such as, pressure, kinetic energy, density and humidity, thus making the model 6.4. Satellite images This section presents a real case, where the performance of the proposed model is compared with the evolution of a clouds from satellite images. The images came from INPE (Instituto Nacional de Pesquisas Espaciais), Brazil. The tests used images from 6 p.m. to 8 p.m. Fig. 11 shows the original images used for the experiment. Fig. 13 shows the cloud images with their respective temperatures in which the scale is given in Fig. 12. Fig. 14 shows two images at 7 p.m.: the result of simulating the proposed model, Fig. 14(a); and the satellite image at 7 p.m., Fig. 14(b). Fig. 14 shows four clouds, Fig. 14(a) shows the simulation of the proposed model from 6 p.m., Fig. 13(a), to 7 p.m., and Fig. 14(b) presents the satellite image at 7 p.m. These four clouds were named top (T), bottom (B), center (C), and right (R). All clouds in Fig. 14(a) tend to homogeneity and disappear because of heat exchange with the atmosphere, since the proposed model uses only one layer of the atmosphere; thus, convection does not act on the clouds. Moreover, the simulated area is isolated. In Fig. 14(a), cloud (T) shows a similar shape to its equivalent one in Fig. 14(b) – the satellite image. Nevertheless, in Fig. 14(b), a cold (blue temperature) cloud appears on the top, different to its equivalent in Fig. 12. Scale of cloud colors. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig. 11. Satellite images used in 4th experiment. 11 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 Fig. 13. Satellite images converted to cloud images. Fig. 14. Result of the first simulation with satellite images. more reliable. Another proposal is to simulate a three-dimensional space, approximating the model of a real system. The proposed model also may be implemented using the parallel computing paradigm, improving its performance and, consequently, its ability to perform in real time. The latter proposal is justified by the increase of variables involved, which feature a real atmospheric system. Thus, parallel computing may increase the model performance in a more complex scenario. Baik, D.-K. (Ed.), 2005. Systems Modeling and Simulation: Theory and Applications. Springer. Dobashi, Y., Kaneda, K., Yamashita, H., Okita, T., Nishita, T., 2000. A simple, efficient method for realistic animation of clouds. In: Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques – SIGGRAPH’00, pp. 19–28. F, F., N, A.-S., E.-S, M., Hassan, A.A., 2010. Modeling and simulation of a single phase grid connected photovoltaic system. WSEAS Trans. Syst. Control 5 (1), 16–25. Fishman, G.S., 1996. Monte Carlo: Concepts, Algorithms, and Applications. Springer Series in Operations Research. Springer-Verlag. Gould, H., Tobochnik, J., 1988. An Introduction to Computer Simulation Methods: Applications to Physical Systems. Addison-Wesley. Guo, H., Pang, J., Shan, Z., 2012. Real-time dynamic cloud simulation based on multicore and multi-thread. In: 2012 Fourth International Conference on Multimedia Information Networking and Security, pp. 241–245. Harris, M.J., 2003. Real-time cloud simulation and rendering. In: ACM SIGGRAPH 2005 Courses – SIGGRAPH ’05, p. 222. Kajiya, J.T., Herzen, B.P.V., 1984. Ray tracing volume densities. In: Proceedings of the 11th Annual Conference on Computer Graphics and Interactive Techniques – SIGGRAPH’84, ACP, pp. 165–174. M, N., P, J.-O., Popescu, M.-C., 2010. Modelling and simulation of traffic lights on road intersections. WSEAS Trans. Syst. Control 5 (7), 529–538. Miyazaki, R., Yoshida, S., Dobashi, Y., Nishita, T., 2001. A method for modeling clouds based on atmospheric fluid dynamics. In: Ninth Pacific Conference on Computer Graphics and Applications, pp. 363–372. Miyazaki, R., Dobashi, Y., Nishita, T., 2002. Simulation of cumuliform clouds based on computational fluid dynamics. In: 2002 Proceedings of Eurographics, pp. 405–410. Nagel, K., Raschke, E., 1992. Self-organizing criticality in cloud formation. Phys. A Stat. Mech. Appl. 182 (4), 519–531. Neyret, F., 1997. Qualitative simulation of convective cloud formation and evolution. In: Animation and Simulation’97, Springer, pp. 113–124. Overby, D., Melek, Z., Keyser, J., 2002. Interactive physically-based cloud simulation. In: 10th Pacific Conference on Computer Graphics and Applications, pp. 469–470. Acknowledgments The authors gratefully acknowledge the Coordination for the Improvement of Higher Education Personnel (CAPES) and Pontifical Catholic University of Minas Gerais (PUC Minas), which supported this work, and the assistance provided by the Center of Climatology of PUC Minas. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.envsoft.2019.104537. References Andrews, D.G., 2000. An Introduction to Atmospheric Physics. Cambridge University Press. 12 A.R. Silva et al. Environmental Modelling and Software 122 (2019) 104537 Pedgen, C.D., Shanon, R.E., Sadowsky, R., 1990. Introduction to Simulation Using SIMAN. McGraw-Hill. Qiu, H., Chen, L.-T., Qiu, G.-P., Yang, H.A.O., 2013. Realistic simulation of 3d cloud. WSEAS Trans. Comput. 12 (8). Rennard, J.-P., 2002. Artificial Life: where Biology Meets Computer Science. FrontMatter & Associates Literary Agency. Rodriguez, R., Lira, J., Rodriguez, I., 2012. Subsidence risk due to groundwater extraction in urban areas using fractal analysis of satellite images. Geofis. Int. 51 (2), 157–167. Rogers, R., Yau, M., 1988. A Short Course in Cloud Physics, third ed. Elsevier, Burlington. Stam, J., 1999. Stable fluids. In: Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques – SIGGRAPH ’99, pp. 121–128. Steiner, J.T., 1973. A three-dimensional model of cumulus cloud development. J. Atmos. Sci. 30 (2), 414–435. Vianello, R.L., 2000. Basic meteorology and aplications. In: Potuguese, Editora UFV, Viçosa. Wolfram, S., 2002. A New Kind of Science. Wolfram Media. Yanagita, T., Kaneko, K., 1997. Modeling and characterization of cloud dynamics. Phys. Rev. Lett. 78 (22), 4297–4300. 13