Monte-Carlo simulations and theoretical modeling are used to study the statistical failure modes and associated lifetime distribution of unidirectional 2D and 3D fiber-matrix composites under constant load. Within the composite the fibers weaken over time and break randomly, and the matrix undergoes linear viscoelastic creep in shear. The statistics of fiber failure are governed by the breakdown model of Coleman (1958a), which embodies a Weibull hazard functional of fiber load history imparting power-law sensitivity to fiber load with exponent ρ, and Weibull lifetime characteristics with shape parameter β. The matrix has a power-law creep compliance in shear with exponent α. Fiber load redistribution at breaks is calculated using a shear-lag mechanics model, which is much more realistic than idealized rules based on equal, global or local load-sharing. The present study is concerned only with the 'avalanche' failure regime discussed by Curtin and Scher (1997) which occurs for sufficiently large ρ, and whereby the composite lifetime distribution follows weakest-link scaling. The present Monte-Carlo failure simulations reveal two distinct failure modes within the avalanche regime: For larger ρ, where fiber failure is very sensitive to load level, the weakest link volume fails in a 'brittle' manner by the gradual growth of a cluster of mostly contiguous fiber breaks, which then abruptly transitions into a catastrophic crack. For smaller ρ, where this load sensitivity is much less, the weakest link volume shows 'tough' behavior, i.e., distributed damage in terms of random fiber failures until the failure of a critical volume and its catastrophic extension. The transition from brittle to tough failure mode for each ρ within the avalanche regime is gradual and depends on the matrix creep exponent α and Weibull exponent β. Also, as α increases above zero the sensitivity of median composite lifetime to load level increasingly deviates from power-law scaling known to occur in the elastic matrix case, α = 0. By probabilistic modeling of the dominant failure modes in each regime we obtain distribution forms and various scalings for damage growth, and for carefully chosen sets of parameter values we analytically extend simulation results on small composites (limited by current computer power) to more realistic sizes. Our analytical strength distributions are applicable for ρ > 2 in 2D, and ρ 4 in 3D. The 2D bound coincides with the avalanche-percolation threshold derived by Curtin and Scher (1997) using entirely different arguments. © 2004 Kluwer Academic Publishers.