  Measuring the 3D distribution of mass on galaxy cluster scales is a crucial
test of the LCDM model, providing constraints on the behaviour of dark matter.
Recent work investigating mass distributions of individual galaxy clusters
(e.g. Abell 1689) using weak and strong gravitational lensing has revealed
potential inconsistencies between the predictions of structure formation models
relating halo mass to concentration and those relationships as measured in
massive clusters. However, such analyses employ simple spherical halo models
while a growing body of work indicates that triaxial 3D halo structure is both
common and important in parameter estimates. The very strong assumptions about
the symmetry of the lensing halo implied with circular or perturbative
elliptical NFW models are not physically motivated and lead to incorrect
parameter estimates with significantly underestimated error bars. We here
introduce a Markov Chain Monte Carlo (MCMC) method to fit fully triaxial models
to weak lensing data that gives parameter and error estimates that fully
incorporate the true uncertainty present in nature. Applying the MCMC triaxial
fitting method to a population of NFW triaxial lenses drawn from the shape
distribution of structure formation simulations, we find that including
triaxiality cannot explain a population of massive, highly concentrated
clusters within the framework of LCDM, but easily explains rare cases of
apparently massive, highly concentrated, very efficient lensing clusters. Our
MCMC triaxial NFW fitting method is easily expandable to include constraints
from additional data types, and its application returns model parameters and
errors that more accurately capture the true (and limited) extent of our
knowledge of the structure of galaxy cluster lenses. (abridged)
