In this work, the application of the cellular automata for the seismic first breaks time estimation is presented. Cellular automata (CA) algorithms are usually described by set of simple rules applied onto the grid of cells, which can represent one of the few different discrete states. In spite of this simplicity, these algorithms can still simulate variety of complex physical processes, i.e. model the fluids or gases behavior, but also simulate and predict forests fire or the spread of diseases (Turcotte 1997). What is more, these algorithms are easy to accelerate using one of the most advanced parallelization techniques, like graphics card programming - the calculation can be then performed even 1000 times faster. In this work one of the popular CA, presented by Hardy et al. (1973), was modified, according to Rothman (1987), to perform the simulation of the seismic wave propagation in geological medium. This CA was initially designed to simulate gas behavior in a reservoir, but it was modified in the following way: each cell of the grid contains set of values, which represents particles - positive and negative - that describes expansion and compression, respectively, of that cell. Then, both particle sets are calculated separately in two steps: the first is the propagation, when each particle is moved to next cell according to its movement; the second is collision, when rules are applied on each particle, if there is more than one particle in the cell. The collision rules can be divided into categories that represent their complexity, i.e. more complex behavior can be achieved if static particles are implemented, but this requires additional rules of collisions with moving particles. All the rules have to preserve the principle of conservation of mass and momentum. The major modification of the algorithm described by Rothman (1987) is the implementation of the Boltzmann lattice method (Huang 2007), which changes the regular rectangular cellular grid for the triangular network. This modification provides more option when implementing algorithm, especially the collisions rules. On the one hand, triangular grid requires slightly different approach to represent it in computer memory, just like the implementations of the rules. On the other hand, triangular grid can produce better simulation results compared to the classical rectangular grid. The results of this simulation can be used in further calculations, i.e. to solve the inverse problem. To make this CA useful, it has to provide either better than similar methods errorless results or significantly shorter computation time with comparable errors. In this work the CA is compared with mathematical wave propagation solution and the Shortest Path Method (SPM). SPM uses graph theory to reconstruct seismic ray trajectory in a real, geological medium. This method was first described in Moser (1991). As it is more time-consuming than linear method, for large model the parallel approach is necessary (Pięta et al. 2013). One of the major advantages is the possibility of optimization of the modeling quality and calculation time by adjusting number of accessory nodes. This method is easily extendable to cover the anisotropy case. This work describes, verifies and validates the CA algorithms to checks the potential of use in further calculations for solving the forward problem in seismic modeling. The research shows that the presented CA method, in spite of its speed and parallelization options, can not be applied if high precision is required.