A four-dimensional non-linear dynamical system resulting from a truncated Fourier representation of the conservation and constitutive equations, for an Oldroyd-B fluid is used to examine the viscoelastic flow in the context of the Rayleigh-Benard thermal convection set-up. Numerical simulations are performed and the results are studied by means of time signature, phase portraits, power spectrum, Poincare map, Lyapunov exponents and bifurcation diagrams. Depending on the flow parameters and initial conditions the asymptotic behavior can be stationary, periodic, quasiperiodic or chaotic. The aim of this paper is to identify the route(s) to chaos, and illustrate the dynamical response of the flow with the change of the control parameters. Fluid elasticity and fluid retardation alter the flow behavior in comparison to inertia-dominated Newtonian flow. Although both flows, Newtonian and viscoelastic, follow similar routes to and from chaos as a function of inertia, elasticity tends to precipitate the onset of chaos. Narrow regions or windows with ordered periodic orbits appear in some regions where the chaotic behavior seems to be dominant. Although the numerical solution can be chaotic, there is a great deal of order, or structure, to deterministic chaos. For example, viscoelastic chaotic attractor is characterized by a higher fractal dimension than the Newtonian chaotic attractor. (C) 2003 Elsevier B.V. All rights reserved.