In proton therapy, cancer patients are irradiated with high energy protons. For a successful treatment it is important that the location with the highest energy deposition, the so-called Bragg-peak, is located inside the tumor and not in the healthy surrounding tissue. Here, we investigate if the iono-acoustic wave field generated by the protons can be used to monitor the Bragg-peak location during treatment. To this end we present a new numerical method to model the pressure field generated by a clinical proton pencil beam. To compute the field, we convolve a 3-D Greens function, representing the impulse response of the medium, with a volume density of injection rate source. This source describes the expansion of the medium due to a local temperature increase caused by the energy deposited by the protons. To image the proton dose distribution, we first measure the pressure field synthetically by a matrix transducer positioned below the Bragg peak in a plane parallel to the beam. Next, we use these measurements to solve the linear inverse problem iteratively. To regularize the inversion, we take the temporal behavior of the dose deposition as prior knowledge. For the presented example, where the pencil beam has a proton range of 63 mm we are able to reconstruct the location of the Bragg peak within 4 mm accuracy.