  We reanalyze gravitational waves from binary-neutron-star mergers GW170817
and GW190425 using a numerical-relativity (NR) calibrated waveform model, the
TF2+_Kyoto model, which includes nonlinear tidal terms. For GW170817, by
imposing a uniform prior on the binary tidal deformability $\tilde{\Lambda}$,
the symmetric $90\%$ credible interval of $\tilde{\Lambda}$ is estimated to be
$481^{+436}_{-359}$ and $402^{+465}_{-279}$ for the case of
$f_\mathrm{max}=1000$ and $2048~\mathrm{Hz}$, respectively, where
$f_\mathrm{max}$ is the maximum frequency in the analysis. We also reanalyze
the event with other waveform models: two post-Newtonian waveform models
(TF2_PNTidal and TF2+_PNTidal), the TF2+_NRTidal model that is another NR
calibrated waveform model, and its upgrade, the TF2+_NRTidalv2 model. While
estimates of parameters other than $\tilde{\Lambda}$ are broadly consistent
among various waveform models, our results indicate that estimates of
$\tilde{\Lambda}$ depend on waveform models. However, the difference is smaller
than the statistical error. For GW190425, we can only obtain little information
on the binary tidal deformability. The systematic difference among the NR
calibrated waveform models will become significant to measure $\tilde{\Lambda}$
as the number of detectors and events increase and sensitivities of detectors
are improved.
