Most drugs have to cross cell membranes to reach their final target. A better understanding of the distribution, interactions, and dynamics of biologically active molecules in model bilayers is of fundamental importance in understanding drug functioning and design. H-2 NMR quadrupole splittings (Delta v(Q)) and longitudinal relaxation times (T-1) from the aromatic ring of benzyl alcohol-d(5) (C-0), a commonly used anesthetic, and a series of linear alkyl benzyl-d(5) ethers with chain lengths from 1 to 12 carbon atoms (C-1-C-12), were measured. The molecules were dissolved in a nematic discotic lyotropic liquid crystal solution made of tetradecyltrimethylammonium chloride (TTAC)/decanol (DeOH)/NaCl/H2O. Values of Delta v(Q) and T-1 from 1,1-dideuteriodecanol (15% enriched) and DHO (H2O with 0.2% D2O) were also measured. Delta v(Q) of DeOH and DHO remained constant throughout the series. The value of Delta v(Q) of the para position of the ring (Delta v(p)) in C, is 30% smaller than the Delta v(p) of C-0. This is attributed to the existence of an H-bond between the alcohol hydroxyl proton and the solvent, which influences the average orientation of the ring. The relaxation data show that T-1o,T-m is always longer than T-1p and both decrease with the increase in alkyl chain length. Molecular dynamics simulations of the experimentally studied systems were performed. The aggregate was represented as a bilayer. The distribution, average orientation, and order parameters of the aromatic ring of the guest molecules in the bilayer were examined. Rotational correlation functions of all the C-D bonds and the OH bond from H2O were evaluated, allowing an estimate of the correlation times and T-1. According to these results all spins relax in extreme narrowing conditions, except DeOH. Experimental and calculated T-1 values differ at most by a factor of 3. However, the order of magnitude and the observed trends are well reproduced by the calculations. The aromatic ring of C-0 possesses a unique average orientation in the bilayer. For the ether series, the orientation is modified and the C-2 symmetry axis of the aromatic ring is exchanging between two orientations averaging the quadrupole splittings from the ortho and meta positions. The simulation supports the existence of an H-bond between C-0 and the solvent not found in the ethers, which should be responsible for the observed differences.