In this paper a numerical method is developed to investigate the tertiary current distribution in an unsteady laminar natural convection multi-ion electrodeposition system. The simulations are based on the dilute-solution theory, with coupling of concentrations of the ionic species through the electric field. The constraint of electroneutrality is applied. For the fluid flow in the electrolytic cell, the conservation equations of mass and momentum must be obeyed. The other governing equations consist of the conservation equation of mass for each ionic species and the conservation equation of electric current. In a finite-volume framework, all the conservation equations are solved by a segregated iterative approach. Also, a general treatment of electrode boundary conditions is proposed and makes it feasible to handle various operational modes of electrodeposition. This capability is important because in many real applications the operational mode is to tune the applied total electric current instead of the applied electrode voltage. The numerical method is validated by an analytical solution from the literature and used to study the copper electrodeposition with a supporting electrolyte in a rectangular cavity of high aspect ratio. Meanwhile, the influence of the operational mode has been studied in detail. For the mode of prescribed total current, results show that the variation of current density along the electrode surface may reach 30% of its average value.